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The procedure for the design of lowpass, highpass 
bandpass. and bandstop digital filters from lowpass 
prototypes is well established, and is based on transferring 
the prototype freguercy response to meet desired 
characteristics . This transformation can be accomplished in 
either the analog (s) domain or the digital (z) domain 
(Figure 1 ) . 
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Figure 1. Two methods for digital filter design. 
Note that each of the horizontal arrows actually 



represent 
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Therefore . 



the two 



classical design techniques. 



analoo- 



dioital and digital-digital design, are based on doing the 
frequency response transformations in their respective analog 
or digital domains. Since the theory for the standard filter 
types (Butterworth, Chebyshev, and elliptic) is developed in 
the analog or 's’ domain, the two methods for designing 
digital filters differ basically in the sequence of the 
design steps. In the analog-digital method the analog 
prototype is transformed and then digitized, while in the 
digital-digital method the analog prototype is digitized and 
then transformed. 

Each of the transformations (arrows in Figure 1) involves 
substituting a function of 's' or ’z’ for the appropriate 
variable in the filter transfer function . The purpose of 
this thesis is to develop efficient algorithms to carry out 
these transformations and then to implement these algorithms 
in a computer program to aid in the design of digital 
filters . 

Computer-aided design of digital filters offers a 
significant advantage since it allows the digital filter 
coefficients to be calculated guicklv and accurately. The 
accuracy of the filter coefficients is an especially 
significant point since the frequency response of the final 
filter is very sensitive to small inaccuracies in the 
calculated coefficients. As an example, compare the 
freouency response for the filters in Figure 2 and 3. Fioure 



2 shows the frequency response magnitude for a 1/2-dB ripple 

% 

Chebyshev bandpass filter, and Figure 3 shows the frequency 
response when one of the coefficients has been changed by 
1/1000 th . - 




Figure 2. Chebyshev Bandpass Filter 

. 00229 (2‘ -1) :! 

H ( r. ! = 

z ( - . 8 7 7 7 f; 7 ' 8 0 0 7 6 r •'* - 1 . 5389t :< +2 . 46277z 2 6749^7+ . 67529 




Figure 3. Same Filter With Denominator Coefficient of 

z changed to -.67594 
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B. GENERAL DISCUSSION 



Each of the transformations involves substitutina a ratio 
of two polynomials 'of s or z into the filter transfer 
function which is also a ratio of two polynomials in s or z . 
namely 

Y ( y ) I 

H ( y ) = j p ( x ) 

U (y ) I y= (1.1) 

a (x) 

The variables 'x' and 'y' can be either 's’ or 'z' depending 
on the transformation to be performed. Therefore, we have 
either 



ao + ai y + a2 v 2 + 

H (y ) = : 

bo + bi y + b2 y 2 + 



+ a>r y m 



+ bn y" 



P ( x ) 

ly= 

g (x) 



( 1 . 2 ) 



or 



ao + ai 



H ( x ) = 



bo + bi 



p ( x ) 



g (x) 



p ( x) 

+ 

g (x) 



(p (x ) ) 2 
a2 + 

(a (x ) ) 2 



(p (x) ) 2 

bi + 

( g ( x ) ) 2 



(d ( x ) ) ® 

+ am 

(a (x) ) m 



(p (x) ) n 
+ bn — — 

( g ( x ) ) " 



( 1 . 3 ) 



In general, the order of the numerator 'm', 
order of the denominator 'n' (except for the 
in the diaital domain where the result 



transformation is to make n=m) 



We can 



is less than the 
transformations 
of the bilinear 
generalize the 
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results, however, by assuming that the order or the numerator 



polynomial Y(x) is equal to the order of the denominator 
polynomial U(x) and making the additional coefficients 

am+i , am ♦ 2 a n equal to zero. Then multiplying both the 

numerator and denominator of the above equation by (g(x))" 
gives 

ao(g(x)) n + a< p (x ) ( g (x ) ) n _ 1 + . . . + a n (p(x)) n 

H (x) = — (1.4) 

bo i g ( x ) ) n + bj p ( x ) ( g ( x ) ) r * ~ 1 + . . . + bn ( p ( x ) ) D 

or, in a different form. 



i=n 

I ai p (x) 1 g (x) D - ’ 
i=0 

H ( x ) = . (1.5) 

i=n 

Z b, p (x) 1 g (x) 11 - 1 
i = 0 



Now with the numerator and denominator of the same form, 
we can use the same transformation for both the numerator and 
denominator polynomials of the filter transfer function, 
namely , 



a’ (x)= an '♦ ai ' x+ a? ' x 2 + . 



. + 



a P ' x n ' 



( 1 . 6 ) 

i=n 

= Z a) p (x) 1 g (x) D * 1 
i = 0 



In Equation 1.6, the a> are the coefficients of the old 
numerator or denominator polynomials and ai ’ are the new. 
transformed numerator or denominator polynomial coefficients. 
Since the new coefficients a* ’ , are a linear combination of 
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the old coefficients ai , and the coefficients of p(x) and 
a(x), this suggests the use of a matrix transformation; 

a ' = aA (n) (1.7) 

where 

a is the coefficient row vector for the old 
numerator or denominator polynomial of 
length=n+l . 

a = [an 3i 32 ... 3n 1 



a' is the coefficient row vector for the new 
numerator or denominator polynomial of 
length=n ' +1 . 



a ' = [ a ’ o a ' i a ' ? . . . a n ’ ] 



A(n) is the transformation matrix with dimensions 
(n+1 ) x ( n ' +1 ) . The elements of A(n) depend upon the 
coefficients of p(x) and o(x). Note that if the 
rows and columns of A(n) are numbered 0 to n ’ , then 

i=n 

a • ' = I a, A (n) , i (1.8) 

i = 0 



(where A(n)i i denotes the element in row i, 
column j of matrix A(n) .) 

n is the order of the transfer function before the 
transformation 

n' is the order of the transfer function after the 
transformation. Note that n’ equals n times the 
highest order of p(x) and g(x). 

For example, when n=l and n ' =2 , p(x) and g(x) are second 

order and Equation 1.7 becomes. 



1 ar 


S • ' 1 A ( 2 ) o n 


A ( 2 ) o i 


A ( 2 ) o 2 i 




i 


i 


— j i 

i A ( 2 ) : f, 


A (2) i , 


j = 

A ( 2 ) « 2 I 

j 


iaO ' 

i 


al ' a2 ' 1 

i 

_) 

(1 
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In the case of the lowpass digital, highpass digital, and 
bilinear transformations, the two functions p(x) and g(x) are 
first order. Professor P.H. Moose of the Naval Postgraduate 
School of Monterey has developed the matrix A(n) for these 
two cases by expanding Equation 1.6 using the binomial 
theorem (Reference 1). This method cannot be used for 
transformations involving higher order polynomials such as 
the bandpass and bandstop transformations and it is difficult 
to develop because intricate manipulation of the binomial 
coefficients are involved. The binomial expansion approach, 
therefore, was not pursued. 

A general method has been found that develops the 
required transformation matrices in an iterative manner. In 
other words, the matrix used to transform a numerator or 
denominator filter polynomial of order n is developed from 
the matrix for the polynomial of order n-1. This approach is 
easy to implement in a computer program and can be used for 
any of the polynomial substitutions needed. The next chapter 
is concerned with developing the transformation matrices for 
each specific transformation. However, as an introduction to 
tne method used, consider the generic case where both p(x) 
and g(x) are second-order polynomials. Since all of the 
analog and digital transformations involve substitutions with 
polynomials of second order or less, this second-order 
Generic case will be all that is needed to develop the 
specific transformations later. 
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With p(x)=ax 2 + bx + c , g(x) 



dx 2 + ex + f , and 



substituting into Equation 1.6 for the first order case of 
n=l (remembering that with p(x) and g(x) of second order, the 
resulting polynomial will be second order) , we get 

a' (x) = ao ’ + ai 'x + a 2 ’ x 2 = an (dx 2 + ex + f ) + at (ax 2 +bx+c) 

( 1 . 10 ) 



Using Equation 1.8 or 1.9 it is easy to pick out the matrix 
coefficients A(1 )j i (i=0,l; j=0,l,2) for the first order 
generic case. The coefficients are, 





i — 




' 






If 


e 


d j 


(1.11) 


Ad) = 


1 c 


b 


a i 






L_ 




_J 







Now with r.= 


2 (n ' 


=4) substituting into 


Equation 1.6 gives, 




a ' ( x ) = ao ' 


+ ai 


' x + a 2 ' x 2 + a3 ’ x 2 + 


a4 ’ x 4 (1.12) 




= ao ( dx‘ +ex+f ) 2 


+ ai (ax 2 +bx+c) (dx 2 +ex+f) +a 2 (ax 2 +bx+c) 2 


and 


the matrix 


A ( 2 ) 


is ; 




A ( 2 ) 


i f 2 
= ifc 

1 c 2 

L_ 


2ef 

ec+bf 

2bc 


2df +e 2 2de 

cd+fa+eb db+ea 

2ac+b 2 2ab 


d 2 i 

dai (1.13) 

a i 

_J 



To see how this matrix is derived from the previous one, 
consider the first row of the matrices in Equations 1.11 and 
1.13. For both matrices, this row is multiplied by the a 0 
coefficient in the old polynomial. For n=l we have a 0 p(x) 
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and we pick out the coefficients of the powers of x (x- 1 ) to 



get the matrix elements A(n) 0 t . For n=2 , we have a y (p(x) ) 2 
and we do the same thina only now we have the additional 
factor of p(x). The matrix elements can be calculated by 
hand in this fashion but it requires a lot of algebra. A 
simpler method is to use the matrix elements of the previous 
matrix and get the appropriate elements by convolvino 
coefficients. For example, consider the sequence of elements 
in row 0 (the top row) of A(l) , namely 



(1.14) 



To get the elements for row 0 of A(2) . we need to multiply 
(convolve coefficients) by another p(x). Therefore, with 
rows numbered C to n, and columns number 0 to n' , we have the 
following . 



For row 


0 column 


0, 


the 


appropriate 


matrix element is 


A ( 2 ) o o = 


d e 


f 






OJ 

II 






f 


e 


d 


(1.15) 


For row 


0 column 


1 , 


the 


appropriate 


matrix element is 


A ( 2 ) o i = 


d 


e 


.L 




= 2fe 






f 


e 


d 


(1.16) 


For row 


0 column 


2 , 


the 


appropriate 


matrix element is 


A ( 2 ) o ? — 


d 


d 


e 


f 


= 2f d+e 2 






f 


e 


d 


(1.17) 
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For row 0 column 3, the appropriate matrix element is 



A ( 2 ) o 3 = 



For row 0 



f 

column 4 . 



d e f = 2de 

e d • (1.18) 

the appropriate matrix element is 



A(2)o<= d e f = d 2 

fed (1.19) 

This convolution of the coefficients of the old matrix with 
the coefficients of p(x) can be carried out for every row of 
the old matrix to obtain all the rows of the new matrix 
except the last row. In the case of the last row, we are 
multiplyina by an additional a(x) rather than an additional 
p(x) and so must obtain the last row of the new matrix by 
convolving the coefficients of the last row of the old matrix 
with the coefficients of g(x). 

In summary we have , 



(i) for row 0 to row n-1: 



-row i of A(n) = row i of A(n-l) * coefficients of p(x) 

( 1 . 20 ) 



(ii) for row n: 

-row n of A(n) = row n-1 of A(n-l) * coefficients of g(x) 

( 1 . 21 ) 

where * indicates linear convolution of the coefficients. 

To see how this method of convolving coefficients works, 
consider multiplying two polynomials: 



1 0 



f 1 ( x ) 



f 2 ( x ) f 3 ( x ) 



( 1 . 22 ) 



Since 'x' is the independent variable, no matter what letter 
we use, the polynomials will remain the same. Substituting 
’ z" 1 ' for 'x' in the above equation gives us, 

f 1 (z- 1 ) = f 2 (z- 1 ) f 3 (z- 5 ) (1.23) 

and we now have the familiar form where we can either 
multiply the two functions of ' z' or convolve the sequences 
found by taking the inverse z-transform. For f2(z -1 ) and 
f 3 ( z - 1 ) we f ind 



f 2 (z- 1 ) 


— clo 


+ ai z - 1 


+ a 2 z - 2 + . 


. . + a n z - n 




f 3 (z- 1 ) 


II 

tr 

o 


+ bi z - 1 


+ b 2 z - 2 + . 


. . + bm Z - ® 


(1.24) 



It is now clear that we can find the coefficients of the 
product fl(x) by convolving the sequence { a n i with the 
sequence 1 bm I to produce the sequence f Ca I where, 

I Bn ( = I 3o , 3l , 32 , 33 , , 3n 1 

! bm I ™ { bo , Id i , Id 2 , b3 , . . . , bm I 

and with q -= m + n (1.25) 

f Cq I = f Co . Cl , C2 , C3 Cq > 

Then by takino the z-transform of the sequence 1 Co I and 

% 

substitutina 'x' for ' z ~ 1 ' we get, 

f 1 ( x ) = Co + C)X + . . . + Cq x^ (1.26) 
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We now have a simple straight-forward method of 
developing the necessary transformation matrices for digital 
filter design by picking out the first-order transformation 
matrix using Equations 1.6 and 1.11. and then using the 
convolution of coefficients method iteratively on the rows of 
each successive transformation matrix with Equations 1.20 and 
1.21 until the appropriate matrix A(n) has been generated. 
The next chapter will develop the transformation matrices for 
each of the following transformations: 
bilinear transformation 



diai 


tal 


lowpass 


prototype 


to 


digital 


lowpass 


f ilt 


er 


digi 


t al 


lowpass 


prototype 


to 


digital 


highpass 


f il 


ter 


digi 


tal 


lowpass 


prototype 


to 


digital 


bandpass 


fil 


ter 


digi 


tal 


lowpass 


prototype 


to 


digital 


bands top 


f il 


ter 



analoo 


lowpass 


prototype 


to 


analog 


lowpass 


filter 


analog 


lowpass 


prototype 


to 


analog 


highpass 


filter 


analoc 


lowpass 


prototype 


to 


analog 


bandpass 


filter 


analog 


lowpass 


prototype 


to 


analog 


bands top 


filter . 


We will see 


that s 


one of the 


se 


transf 


ormat ions 


are similar 


and consequently. 


we wi 


11 


not 


have 


to develop 


separate tra 


nsf orr.a 


tier, matric 


es 


for ea 


ch case . 
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II . 



DEVELOPMENT OF SPECIFIC TRANSFORMATIONS 
A. INTRODUCTION 

The individual transformations are discussed in the 
following sections. The specific form of the transformations 
used in each case are from Reference 2. Since the final 
result will be to develop a computer program to assist in the 
calculation of digital filter coefficients, confining the 
specific transformations to the form in Reference 2, will 
allow the computer program to be used in conjunction with 
Reference 2 as an aid in designing digital filters. 



B. BILINEAR TRANSFORMATION 



The bilinear 


transformation may 


be used 


to transform an 


analog transfer 


f unct i on . H ( s ) , 


to a 


digital transfer 


function H ( z ) . 


The transformation 


from 


H ( s ) to H ( z ) is 


described by 








H ( z ) = 


H(z) 1 








i z — 1 




(2.1) 




1 s = 








z + 1 






Referring to Equ 


ation 1.1, p ( x ) = z - 1 


and g(x) 


= z + i. Since both 



of the transformation polynomials p(x) and g(x) are first 
order, the bilinear transformation will not change the order 
of the transfer function. Therefore, for an analog transfer 
function of order 'n' , the transformation matrix A(n) will be 
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(n+1) by (n+1) . From Equation 1.6, the numerator or 

denominator polynomial of the transfer function H(z) is given 
by 

i=n 

a ’ (z ) = ao 1 + a, ' z+ . . . + a n ’ z n = I ai ( z-1 ) 1 ( z + 1 ) n - 1 

i=0 (2.2) 



where the unprimed coefficients, a* , come from the analog 
transfer function’s numerator or denominator polynomial. In 
matrix form, with n=l and referring to Equation 1.11, the 
first order matrix A(l) for the bilinear transformation is 



A ( 1 ) 




-1 



i i 
1 1 



(2.3) 



Usina the 
Chapter 1 , 
from the 
tr ans f orma t 



convolution of 


coefficients 


method de 


the 


transformation 


matrix A(n+1) can be 


ma 


trix A ( n ) . 


The results 


for the 


ion 


matrices A(l) 


through A(5) 


are give 



veloped in 
developed 
bilinear 
n in Table 



1 . 
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FIRST-ORDER 



ACL) = 



A ( 2 ) = 



A ( 3 ) = 



A ( 4 ) = 



A ( 5 ) = 



TABLE 1 

TO FIFTH-ORDER BILINEAR TRANSFORMATION MATRICES 



1 I 

_J 



1 


2 


1 i 

I 








-1 


0 


1 1 
I 








1 


-2 


1 

1 1 








L_ 




_J 








l 


3 


3 


1 

i 


i 

i 




-l 


-1 


1 


1 


i 




l 


-1 


-1 


1 






-l 


3 


-3 


1 






— 






_l 






H i 

i 


4 


6 


4 


1 ! 

i 




i 

|-i 


-2 


0 


2 


i 

1 ! 

i 




i 


o 


-2 


0 


i 

i 1 

i 




-i 


2 


0 


-2 


i i 
1 




i 


-4 


6 


-4 


i 

1 i 




— 








_! 




— 










— 1 


i 


5 


10 


10 


5 


1 


-i 

i 


-3 


-2 


2 


3 


1 

« 


i i 
1 


1 


-2 


-2 


1 


1 


i 

i-i 

1 


1 


2 


-2 


-1 


1 


i i 

i 


-3 


2 


2 


-3 


1 


i 

i-i 


5 


-10 


10 


-5 


1 


!_ 










_l 






Referring to Table 1 and keeping in mind the convolution 
of coefficients method that produced the matrices, the 
following relations are apparent and will be useful in 
computer implementation. 

A(n)o» (the first row) = ( ” ) for j=0 . 1 . . . . , n (2.4) 

where ( j ) indicates the binomial coefficient 

A(n) i o (the first column) = (-1)‘ for i=0,l,...,n (2.5) 

A(n)tn (the last column) = 1 for i=0,l,...,n (2.6) 

A(n)<n-i>t = (-1) n " J A (n) i t for i=0 , 1 ..., integer (n/2 ) +1 (2.7) 

A(n)ii = A ( n-1 ) j j + A ( n-1 ) i < i - j » for i=0.1,...,n-l (2.8) 

j=l , 2 , . . . , n-1 

Notice that the values of ' the matrix elements for the 
bilinear transformation are whole numbers. This is important 
from, a numerical accuracy standpoint since whole numbers can 
be expressed precisely in floating point form. This means 
that the bilinear transformation will cause a loss of 
accuracy of the original coefficients only because of the n+1 
multiplications and additions used to calculate each new 
digital coefficient. 

C. LOWPASS TO LOWPASS DIGITAL TRANSFORMATION 

The lowpass to lowpass diaital transformation is used to 
transform a lowpass prototype digital filter H(z) into a 
lowpass digital filter H(z) . The polynomial substitution 
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used is 



H(z)= H P (z) 



z-a 

z= 

1-az 



(2.9) 



where the subscript p indicates the prototype digital filter 
transfer function. The design constant a is determined 
[References 2 and 3] by 

sin (6r / 2 - 6c ' / 2 ) 

a = 

sin ( ©c / 2 - ©r " / 2 ) (2.10) 
with ©r = the prototype critical frequency (usually pi/2) and 
©r '= the desired critical frequency of the lowpass filter. 
Again, referring to Equation 1.1 with p(x)= z-a and g(z)= 1- 
az , the lowpass digital transformation matrix will not change 
the order of the transfer function and the transformation 
matrix A(n) will be (n+1) by (n+1). Using Equations 1.6 and 
1.11 as before, the first order transformation matrix A(l) 
for the lowpass digital transformation is: 



A (1) 




-a 

1 



_J 



( 2 . 11 ) 



Once again using the convolution of coefficients method, the 
higher-order transformation matrixes A(n) are developed from 
the lower-order matrix A(n-l). The results are shown in 
Table 2 for the first-order to fifth-order transformation 
matrices . 
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A ( 1 ) 



A ( 2 ) 



A ( 3 ) = 



A ' 4 ' = 



TABLE 2 

FIRST-ORDER TO FIFTH-ORDER LOWPASS DIGITAL 
TRANSFORMATION MATRICES 



1 1 
1 


-a 












b 


n 












! 1 
i 


-2a a 2 i 

1 








j -a 
j 


( 1 + a 2 ) -a 


• 








I 

i a 2 

i_ 


1 

-2a 1 i 

_J 








1 


-3a 


3a 2 -a 3 


1 

1 






1 -a 


( l+2a 2 ) ( 


-2a-a 3 ) a 2 


1 

I 






ia 2 


( -2a-a 3 ) ( l + 2a 2 ) -a 


i 

I 

1 






1 

l-C 3 

j 


3a 2 


-3a 1 


♦ 

1 

1 

i 

_J 






1 2 
i 




— 4 a 


6 a 2 


— A ry 3 


a 4 1 

j 




1 -a 

i 


( 1 


*3a 2 ' 


(-3a- 3a 3 > 


( 3a 2 +a 4 ) 


-a 3 1 
| 




1 -v 2 

■ j . 

| 


i ~ 

\ 


2a-2a 3 


• ( 1 4 a 2 -‘•a 4 


^ f -2a-^a 3 ) 


1 

r* 2 I 

w ‘ t 

1 




( — rv 3 
1 


( 

<> w 


c 2 ) 


( — 3 a — 3 a. 3 ^ 


( l + 3a 2 > 


I 

-a { 

] 

1 




1 rv A 

i_ 




-4a 3 


6 a 2 


-4a 


I 

i j 

_i 




! 1 


- 


5a 


10 a 2 


-10a 3 


5a 4 


-a 5 | 
j 


1 -n 

1 


( 1 + 3a 2 > 


(-4a-6a 3 ) 


( 6a. 2 +a 4 ) 


( -4a 3 -a 3 ) 


a4 i 

j 


1 

1 a 2 
j 


(-2a 


-3a 3 ) 


( l + 6a 2 + 3a 4 ) 


( -3a-6a 3 -a 5 ) 


(3a 2 + 2a« ) 


-a 3 j 

j 


1 -or 

1 


( 3a 2 +2a« ) 


( -3a-6a 2 -a 5 ) 


( l + 6a 2 + 3a 4 ) 


( -2a-3a 3 ) 


1 

a2 1 
I 


1 a 4 
i 


(-4a 3 -a 5 ) 


(6a 2 + 4a 4 ) 


( -4a-6a 3 ) 


( l + 4a 2 ) 


-a | 
i 


1 

j ~a.' J 


5a 4 


-10a 3 


10a 2 


-5a 


i i 



l_ —I 
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Referring to Table 2, the following relations are 
observed and will be used in the computer implementation. 



A (n) o j 


(the first 


row) = 


«; 


) (-a) J for j=0,l, ,n 


(2.12) 


A ( n ) * o 


(the first 


column ) 


= 


(-a) 1 for i-0 , 1 , . . . ,n 


(2.13) 


A (n) i j 


= A ( n ) < n “i ) 


( n - ) ) 


for 


i , j=o ,1 n 


(2.14) 


A ( n ) i i 


= A (n-1 ) ( 1 - 


1 ) ( ) - 1 ) 


+ 


( -a) A (n-1 ) ( i - i > i 


(2.15) 



for i=l , 2 . , n-1 

j = l , 2 n 



Note that the elements of the lowpass digital 
transformation matrices are not whole numbers as in the 
bilinear transformation case, but are polynomials of the 
design constant a. Numerical error is introduced in this 
case by both the n+1 additions or multiplications, and the 
error in calculating the matrix elements that multiply the 
original coefficients. 



D. LOWPASS TO HIGHPASS DIGITAL TRANSFORMATION 

The lowpass to highpass digital transformation is used to 
transform a- lowpass prototype to a highpass filter. The 
polynomial substitution for this transformation is: 



H ( z ) 



H P (z) 



z 



z-a 



1-az 



(2.16) 



Since this transformation is the same as that for the lowpass 
to lowpass digital transformation except for the minus sign. 
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we can make the substitution of ’-z' for ’z' in the original 
prototype transfer ' function and then use the same 
transformation matrices as the lowpass to lowpass case. For 
the lowpass to highpass transformation, the design constant a 
is (References 2 ana 3) 



cos (©c / 2 - ©e ’ / 2 ) 

a = (2.17) 

cos (0r /2 + ©c ' / 2 ) 



with ©c - the prototype critical frequency (usually pi/2) and 
0c' - the desired critical frequency of the highpass filter. 



E. LOWPASS TO BANDSTOP DIGITAL TRANSFORMATION 



The polynomial substitution which transforms a digital 
lowpass prototype to a bandstop digital filter is (References 
^ and o ) t 



2a 1-k 

H(z) = Hr (z) I z 2 - z + 

I 1+k 1+k 

I z = 

2a 1-k 

1 - z + z 2 

1+k 1+k 



with 



cos (G,. 1 /2 + ©i * / 2 ) 

a = 

cos (©u ' 12 - ©i ’ / 2) 



(2.18) 



(2.19) 



ana 



K = tan(©r /2)tan(0 u ’/2-0i 72) 



( 2 . 20 ) 



with once again, © c = the prototype critical frequency 
(usually pi/2) , ©u ’ = the desired upper critical digital 
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frequency, and ©i ’ 



the desired low critical diaital 



frequency . 

By letting b= 2a/(k+l) and c= ( 1-k ) / ( 1+k ) . Equation 2.18 can 
be simplified to 



H (z) 



H P ( 2 ) 



2 



2 2 “ b2 + C 



1 - b2 + C2 2 



( 2 . 21 ) 



Again, using Equations 1.6 and 1.11 the first-order 
transformation matrix is 

A(l) = il -b 
i 

i 

j c -b 

L 

Mote that since the substitution polynomials are second 
order, this transformation will double the order of the 

matrix A(n) will be (n+1) by (2n+l). The higher order 
transformation matricies are again developed recursively and 
the results for the first-order to third-order 
transformations are given in Table 3. 



( 2 . 22 ) 
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TABLE 3 

FIRST-ORDER TO THIRD-ORDER BANDSTOP 







DIGITAL TRANSFORMATION 


MATRICES 


A ( 1 ) = 


l 


-b c 










c 


-b 1 








A ( 2 ) = 


1 


-2b 


( b 2 + 2c ) -2bc 


c 2 i 






c 


(-b-bc) 


( 1+b 2 +c 2 ) (-b-bc) 


c 






c 2 


-2bc 


(b 2 +2c ) -2c 


1 




A (3) = 


1 


-3c 


(3b 2 +3c ) 


(- 


-b 3 -6bc ) 


i 


c 


( -2bc-b ) 


( l+2b 2 +bc +2c> 


(-b 3 - 


-2b-2bc-2bc 2 ) ... 

1 


1 

i 

i 


c 2 


(-2bc-bc 2 ) (c 3 +2b 2 c+2c+b 2 ) 




• • • 


i c 3 

L_ 


-3c 2 b 


(3b 2 c + 3c 2 ) 




• » • 



(Redundant entries not shown for A(3) due to space 
limitations. Entries not shown can be found using Equation 
2.25) 

Notice the increased amount of calculation needed to 
determine the matrix elements for the bandstop case. Each 
matrix element is a function of 'b* and 'c' which are 
themselves functions of the design constants a and k. Table 
3 only gives the transformation matrices for first- to third- 
order. However. to illustrate the amount of calculation 
required, the matrix element A(5)25 is 

A ( 5 ) 2 rv = ( -12bc 2 -6bc 3 -8b 3 c-6bc-6b 3 -3b-3bc 4 -6b 3 c 2 -b 5 ) (2.23) 

Since each coefficient in the final filter is found by 
TT'Ul tiplving n + 1 matrix elements by the n+1 original 
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coefficients and then addina the n + 1 products, the accuracy 
of the design constants and the original coefficients can be 
critical in the calculation of the final filter coefficients. 
The symmetry relations for the bandstop transformation 

are: 



A(n)i o (the first column) = c 1 for i=0,l,...,n 



A(n)i.i = A(n)<n-i>( 2 n-i> for i = 0 , 1 n 

j = 0 .1 2n 



A ( n ) i j = A ( n-1 ) i i - bA (n— 1 ) i < j - i > + cA< n-1 ) i < i - 2 > 

for i = 0 ,1 n-1 

j=0 , 1 , . . . , 2n 



(2.24) 

(2.25) 



(2.26) 



F. LOWPASS TO BANDPASS DIGITAL TRANSFORMATION 

Again, the polynomial substitution which transforms the 
lowpass prototype to a bandpass filter is 



H ( z ) = H P ( z ) 



2ak k-1 

z 2 - z + 

k+1 k+1 



2ak k-1 

1 - z + z 2 

k+1 k+1 



(2.27) 



where in this case k = tan (0/2 ) cot (©,. ’ /2 - ©1 ’ / 2 ) and a is 
found from the same expression as the bandstop case. With 
the simplification of b=2ak/(k+l) and c= ( k-1 ) / ( k+1 ) , Equation 
2 . 27 becomes 
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H (z) 



H P (z) 



z 2 - bz + c 

z = - (2.28) 

1 - bz + cz 2 



Since this is the same form as the bandstop 
transformation, we can use the matrices previously developed 
if '-z' is first substituted for ' z' as was done in the 
lowpass and highpass situation. 



G. LOWPASS TO LOWPASS ANALOG TRANSFORMATION 

The lowpass to lowpass analog transformation is used to 

transform an analog lowpass prototype to a lowpass analog 
filter. The polynomial substitution is a simple one, namely 



H ( s ) = Hp (s) 



s = s/w ( 

i 



(2.29) 



where w r is the prewarped critical frequency (w r =tan (0c- /2 ) 
and ©c is the desired digital critical frequency 
( ©r =2 ( pi ) f c / f s ) . This transformation matrix is found to be 



II 

i 

i o 
i 

A ( r. ) = j C 

I . 

i- 

i . 
i 

I o 

L_ 



o o 

1 /Wr 0 

0 1 / Wc 2 

0 0 



0 

0 

0 



1/Wr 0 I 




or A(n) = diaa( 1, l/w c ,1/wo 2. 
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, l/w c n) throughout. 

(2.30) 



Note that in this case, 
depends only on multiplying 
(Wr ) _1 . Therefore, any loss of 
transformation is due to this 
error introduced in calculation 



each new filter coefficient 
the old coefficient of s' by 
accuracy due to the 
one multiplication and the 
of the matrix element (w c ) _1 . 



H. LOWPASS TO HIGHPASS ANALOG TRANSFORMATION 

The polynomial substitution for the lowpass to highpass 
transformation is 



H ( s ) = H P (s) 



S = Wr /s 



(2.31) 



Again using Equations 1.6 and 1.11, the first-order 
transformation matrix is 



1 I 

(2.32) 

0 I 

_j 

The method of convolution of coefficients can be used to get 
the hioher-order transformation matrices. However, notice 
that Equation 2.32 is the same matrix as the lowpass to 
lowpass transformation with the reciprocal taken of each non 
zero element, and then each row reversed. This can be shown 
to be the case in general , so that the lowpass to highpass 
transformation matrix is given by 



A ( 1 ) = 



0 

w c 
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0 



0 



(2.33) 



0 

1 0 0 

I . 

A ( n ) = i . 



w. 



1 I 

i 

r\ I 



0 Wr n - 1 

Wr" 0 



0 0 i 

i 

o o i 



Again. each new filter coefficient depends only on one 
multiplication although the calculation of the matrix element 
itself will require calculations involving raising the design 
constant, w c , to the (n-l) th power where i is the power of s 
in the final filter transfer function. 



I. LOWPASS TO BANDPASS ANALOG TRANSFORMATION 

The polynomial substitution for transforming a lowpass 
analog prototype to a bandpass analog filter is 



H ( s ) = H P (s) | 

I s 2 + wo 2 

! s = 

Bs 



(2.34) 



where wo is defined as the geometric mean of the passband; 
wo = (wuwi ) s , and B is the width of the passband, B=w u -wi . 

(w» and wi are the upper and lower critical frequencies 
respectively). Again using Equations 1.6 and 1.11 and 
picking out the appropriate matrix elements, the first-order 
transformation matrix is given by 
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Ad) 



0 B 0 I 

I 

Wo 2 0 1 I 

_! 



(2.35) 



The convolution of coefficients technique is again used to 
generate the higher order transformation matrices with the 
results for the first-order to third-order matrices listed in 
Table 4. As with the bandpass digital transformation, the 
analoc bandpass transformation will double the order of the 
original transfer function since the substitution polynomial 
is of second order. 



TABLE 4 

FIRST-ORDER TO THIRD-ORDER BANDPASS 
ANALOG TRANSFORMATION MATRICES 



i 0 B 

Ad ) = I 

I wo 2 0 

i_ 



0 ! 

i i 



! 0 



o 



B 2 0 0 j 



A ( 2 ) = iO 

I 



Bwo 2 0 



B 



0 ! 



Wo 4 0 

I 

I 



2wo 2 0 1 i 

J 



A ( 3 ) 




Ud 6 



I 



0 0 

0 B 2 Wo 2 
Bwo 4 0 

0 3Wr 4 



B 3 0 

0 B 2 

2Bwo 2 0 

0 3w~ 2 



0 

0 

B 

0 



i 

o i 

o i 
i 

i! 



Note that each matrix is slightly over half empty and 
that the nonzero elements are simple functions of the desian 
constants Wo and B. 
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J. LOWPASS TO BANDSTOP ANALOG TRANSFORMATION 

For the analog bandstop transformation, the polynomial 
substitution is the reciprocal of the bandpass case: 



H (s) = Ho (s) 



Bs 



( 2 . 36 ) 



s - 



S 2 + Wo 2 



Usina the same method as before, the first-order 
transformation matrix is found to be: 



| Wo 



A ( 1 ) = I 



I 0 

L_ 



0 

B 



il 



( 2 . 37 ) 



Note that this is the same matrix as the bandpass case but 
with the rows in reverse order. Using the convolution of 
coefficients technique the higher-order matrices are 
developed and the results for the first-order to third-order 
transformation matrices are shown in Table 5. Notice that the 
reversal of the order of the rows compared with the bandpass 
case is true in general. Therefore, once again the matrices 
are over half empty and nonzero elements are simple functions 
of B and Wo . 
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TABLE 5 

FIRST-ORDER TO THIRD-ORDER BANDSTOP 
ANALOG TRANSFORMATION MATRICES 



A ( 1 ) 




0 

B 



1 i 

I 

C I 



Wo 4 0 2wo 2 



0 



1 



A ( 2 ) 



0 



Bw 0 4 0 



0 



0 B 2 



B 

0 



0 



0 I 



Wo 6 0 2wo 4 



0 3wo 2 2 0 



1 



0 



A ( 3 ) = 

i 0 



Bwo 4 0 2Bw 0 2 0 

0 B 2 wo 2 0 B 2 



B 0 

0 0, 



i 0 0 0 

L_ 



B 3 0 



0 



0 



K. GENERAL OBSERVATIONS 

Now that the various transformation matrices for digital 
filter design have been developed, it is possible to compare 
the two methods of diaital filter design mentioned earlier 
and depicted in Figure 1. For this comparison, consider an 
analoc lowpass prototype filter transfer function with 
numerator of order mi, and denominator of order n (m less than 
or equal ton). 

First consider the dioital to digital direct design 
method. Transforming the analog prototype to a digital 
prototvoe usino the bilinear transformation will require m+1 
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multiplications and additions per digital prototype 
coefficient for the numerator polynomial, and n+1 
multiplications and additions per digital prototype 
coefficient for the denominator polynomial. Then, each of 
the digital to digital transformations will require 2 (n+1) 
multiplications and additions since a result of the bilinear 
transformation is to make n = m. Therefore, going from an 
analog prototype to a final digital filter will require 
(3n+m+4) multiplications and additions. If lowpass digital 
prototypes are developed prior to the design process, this 
can be reduced to 2 (n+1) since then only the digital to 
digital transformations will be used in designing the final 
filter . 



Now consider the analog design method for the same case 
of an m 1 h -order numerator and n’ b -order denominator lowpass 
prototype transfer function. For the lowpass and highpass 
case, each new coefficient will require only 1 multiplication 
and addition. For the bandpass and bandstop transformations 
a simple expression for the number of calculations per 
coefficient is difficult to find, however referring to Tables 
4 and 5. the number of multiplications and additions per 
coefficient is always less than half the order of the 
polynomial to be transformed, on the average. This gives 
approximately (n-t-m) /2 for the bandpass to bandstop analoa 
transformations. Finally the analog filter must be converted 
to a dicital filter using the bilinear transformation which 
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will require (n+m+2) calculations per coefficient for the 

highpass and lowpass case and 2(n+m+l) for the bandpass and 
bandstop case. The total nur.be r of calculations in using the 
analog design method then is (n+rn+4) for the high/low pass 
and ( 2 . 5 (n+ir.) +2 ) for the bandpass/stop. These results are 
summarized in Table 6. 



TABLE 6 

NUMBER OF MULTIPLICATIONS AND ADDITIONS REQUIRED 
PER COEFFICIENT TO FIND THE DIGITAL FILTER 
COEFFICIENTS FROM A LOWPASS ANALOG PROTOTYPE 
(WITH AN M TH ORDER NUMERATOR AND N T H ORDER DENOMINATOR) 



Analog Desig n 



high/low pass bandpass/stop 



analoc transformations 


2 

l_ 


I 

| . 5 ( n+m) 


bilinear transformation 


i 

! n+m+2 

! 


1 2 ( n+m) +2 1 

! ..... . ! 


Total 


I 

i n+m+4 

j 


! j 

2 . 5 (n+m ) +2 j 




Digital 
high/low pass 


Design 

bandpass/stop 


bilinear transformation 


| 

n+m+2 


n+m+2m j 

j 


digital transformations 

j 


2 (n + 1 ) 

i 


i 

2 (n+1) i 

i 


1 

Total 

1 


3 (n + 1 ) +m+ 1 | 

! i 


! 

! 

3 ( n+1 ) +m+l I 

j 



31 



Referring to Table 6, the advantaae of one method over 
the other depends upon what type of filter is being designed 
and the order of the prototype transfer functions numerator 
and denominator. Consider designing a fifth-order 
low/highpass Chebyshev or Butterworth filter. In this case, 
the order of the numerator polynomial is zero and the analog 
method will have a slight advantage of 9 multiplications and 
additions versus 12 for the digital design with precalculated 
prototypes (and 19 for digital design without precalculated 
prototypes). However, if a bandpass /stop elliptic filter was 
being designed with n=5, then m=4 and the direct digital 
desian method (with precalculated prototypes) has an 
advantage of 12 multiplications and additions compared to 
approximately 24 for the analog method. 

A final consideration in favor of analog design is the 
complexity of the elements of the transformation matrices. 
Since the elements of the digital transformation matrices 
require much more calculation than the elements of the analog 
transformation matrices or the bilinear transformation 
(compare Tables 1-5), the problem of roundoff error is much 
more serious in the case of the digital to digital design 
method . 

L. IMPLEMENTATION 

The final step is to implement the algorithms developed 
previously in a computer program to assist in the teaching of 
dicital filter design. Since the program is intended to be 



32 



used for instructional purposes, the numerical considerations 
discussed earlier are not considered as important as 
conceptual simplicity. The direct design method was chosen 
for implementation as the design constants can be found 
easily without reference to " prewarping The completed 
program. called ' DFCADD ' (Digital Filter Computer Aided 
Direct Design) , is contained in Appendix A with its own 
documentation. The program follows the procedures set forth 
in Chapter 10 of Reference 2. Prototype coefficients are 
stored within the program for first-order to fifth-order 
Butterworth and Chebyshev (l/2,l,&2dB ripples) filters. 
While second-order to sixth-order elliptic lowpass prototype 
coefficients are calculated within the program for any 
combination of 1/2, 1,2, and 3 dB passband ripple and -20,- 
30 ,- 40 .- 50 ,- 60 . and -70 dB stopband attenuation, with the 
help of a program, developed by Professor D.E. Kirk of the 
Naval Postgraduate School and References 4 and 5. As a final 
option, the user may enter his/her own prototype coefficients 
for up to a tenth-order analog prototype. After calculating 
the digital filter coefficients an option exists for the 
rroorar to create a data file, titled “filter*', with the 



digital 


filter macnitude 


and phase 


as a function 


of 


the 


digital 


frequency from z 


ero to pi. 


This data file 


can 


then 


be used 


with any plottinc 


routine to 


generate ' a plot 


of 


the 


freouenc 


v response of the 


filter . 
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APPENDIX 

DFCADD FORTRAN PROGRAM 



u 



FILE* DFCADD FORTRAN A 



9 

10 



XXXX X X X X X X X X XXXX X X X XX X X XXX X X X X XX XX X XXX X ***** XXX X X X X XXX XXX XXXX * X X 

x THIS PROGRAM ASSISTS IN THE CALCULATION OF DIGITAL FILTER x 
x COEFFICIENTS ACCORDING TO THE PROCEDURE SET FORTH IN CHAPTER * 
* 10 (SUMMARIZED ON PG 695) OF - FIRST PRINCIPLES OF DISCRETE * 
x SYSTEMS AND DIGITAL SIGNAL PROCESSING - BY R.D. STRUM AND X 
X D.E. KIRK (REFERENCED) ) X 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
INTEGER ORDER, TYPE, PASS , 0RDER2 , END, AFORDF 
REAL KPR, LAM, OMEGA( 20) ,V(20),AO(20), BO(20),B1(20) 

COMPLEX Z, NEWZ, ZN, ZD, HOFZN, HOFZD, HOFZ 
DOUBLE PRECISION T EMP , SUMR 

DOUBLE PRECISION R , SO , A01 , A02 , BO 1 , B02 , B1 1 , B1 2 , HO 

DOUBLE PRECISION AO 3 , AO A , B03 , BO A , B1 3 , B1 A 

DOUBLE PRECISION S FREQ , AFC, DFC , AL PHA , PI , PI DA , DFCD2 

DOUBLE PRECISION ALC,AUC,DLC, DUC, RK, B, C, DLCD2 , DUCD2 

DOUBLE PRECISION APROT , BI L I NR, DPROT , L PHPTR , BPBSTR , DFL TR 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

FORMAT OF ARRAYS* 

-APROT CONTAINS THE ANALOG PROTOTYPE 
APROT(NORD,PWR) 

-BILINR IS THE ANALOG-DIGITAL BILINEAR TRANSFORMATION MATRIX 
BILINR(ORDER,ROW,COL) 

-DPROT IS THE DIGITAL PROTOTYPE FROM APROT AND BILINR 
DPROT(NORD,PWR) 

-LPHPTR IS THE LPSHP DIGI TAL -DIGITAL TRANFORMATION MATRIX 
LPHPTR(ORDER,ROW,COL) 

-DFLTR IS THE FINAL DIGITAL FILTER 
DFLTR ( NORD, PWR) 

-BPBSTR IS THE BPSBS DIGITAL-DIGITAL TRANSFORMATION MATRIX 
BPBSTR ( ORDER, RON, COL ) 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
DIMENSION APROT (2,0*10) 

DIMENSION BILINRdO, 0*10,0*10) 

DIMENSION DPROT (2,0:10) 

DIMENSION LPHPTRdO, 0:10, 0*10) 

DIMENSION BPBSTRdO, 0*10, 0:20) 

DIMENSION DFL TR ( 2 , 0 : 20) 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
OPEN( UNI T = 8, FI LE= f FILTER DATA ' , STATUS = 1 UNKNOWN * ) 

XXXXXXXXXXXXXXXXXXXXXX*X***xXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

PI = 3 . 14159265358 98D+00 
SPI =3 .14159 
PID4=PI/4.0D+00 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

CLEAR APROT PRIOR TO LOADING 
DO 1 K= 1 , 2 
DO 2 L=0,10 
APR0T(K,L)=0. 0D+00 
CONTINUE 
CONTINUE 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 



XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

BEGIN INTERACTION 
CALL EXCMS (’CLRSCRN’) 

MRI TE( 6 , x ) ’THIS IS A PROGRAM TO CALCULATE DIGITAL FILTERS USING 1 
NR I T E ( 6 , x ) » THE DIRECT DESIGN METHOD.' 

HRI TE( 6 , x ) 

MRITE( 6 ,x) ’xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx* 

NR I TE( 6 , x ) ' XHARNING-THIS PROGRAM IS NOT USER FRIENDLY*!! x» 
WRITE( 6, X ) »x YOU MUST ENTER VALUES OF THE APPROPRIATE TYPEx* 
NRITE( 6 , x ) ’x (REAL OR INTEGER), AND PROPER RANGE X* 

MRITE( 6 , x) ' XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 1 
WRI TE( 6 , x) 

MRITE(6 , x) * xx ENTER 1 TO CONTINUE xx ' 

READ( 5 , x ) LOOK 
CALL EXCMS ( * CLRSCRN * ) 

NRITE( 6 , x ) ’WHAT TYPE OF FILTER DO YOU WANT?' 

MRITE( 6 , x) ' (1-4 ALLON ONLY UP TO 5TH ORDER PROTOTYPES)' 

WRITE( 6 , x ) ' 1-BUTTERNORTH ' 

NRI T E( 6 , X ) ' 2-CHEBYSHEV WITH 1/2 DB RIPPLE ' 

WRI TE( 6 , x ) » 3-CHEBYSHEV WITH 1 DB RIPPLE ' 

WRITE( 6 , x ) » 4-CHEBYSHEV WITH 2 DB RIPPLE' 



n 



NEW00010 
NEW00020 
NEW00030 
NEW00040 
NEW00050 
NEW00060 
NEW00070 
NEW00080 
NEW00090 
NEW00100 
NEW001 1 0 
NEW001 20 
NEW00130 
NEW001 4 0 
NEW00150 
NEW00160 
NEW00170 
N ENO 018 0 
NEW00190 
NEW00200 
NEW0021 0 
NEW00220 
NEW00230 
NEW00240 
N EW 0 0250 
NEW00260 
NEW 00 270 
NEW00280 
NEW00290 
NEW00300 
NEW00310 
N EW00320 
NEN00330 
NEW00340 
NEW 00 350 
NEW00360 
N ENO 0370 
NEMO 038 0 
N EM 0 039 0 
NEN004 0 0 
NEW0041 0 
N EW 00420 
NEMO 0430 
NEMO 0440 
N EW 00450 
NEMO 0 46 0 
NEMO 04 7 0 
NEMO 048 0 
NEMO 0490 
NEN00500 
NEMO 051 0 
NEMO 0520 
NEMO 0530 
NEMO 0540 
NEW00550 
NEW 00 560 
NEW00570 
NEMO 058 0 
NEW00590 
NEW00600 
NEMO 06 1 0 
NEMO 06 20 
NEMO 06 30 
NEW00640 
NEW00650 
NEW0 066 0 
NEW0067 0 
NEW0068 0 
NEMO 06 9 0 
NEW00700 
NEW00710 
NEW00720 
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C 

14 

15 



C 



C 



WRITEC6,x) 
WRITEC6,X) 
WRITEC6,X) 
WRITEC6,X) 
WRITEC 6 , x) 
WRITEC 6 , x ) 
MRI T EC 6 , X ) 
WRITEC6,X) 

WRITEC 6 , X) 



5- ELLIPTIC CALLOWS UP TO A 6TH ORDER PROTOTYPE) 1 
-PASSBAND RIPPLE 0.5, 1.0, 2.0, OR 5.0 DB' 
-STOPBAND ATTENTUATI ON -20,-50,-40,-50,-60, • 

OR -70 DB ■ 

6- OTHER - ALLOWS YOU TO USE UP TO A 10TH ORDER 1 

PROTOTYPE, HOWEVER YOU MUST ENTER THE ANALOG 1 
PROTOTYPE COEFFICIENTS ' 

XX INTEGER 1-6 XX* 



MUST BE AN INTEGER IN THE RANGE 1-5 UNLESS YOU* 
CHOSE TO ENTER YOUR OWN PROTOTYPE C THEN 1-10)' 

OR YOU CHOSE AN ELLIPTIC FILTER C THEN 1-6)' 

OR ENTER "0" IF YOU WANT HELP SELECTING THE ORDER' 
BY CHANGING SPECIFICATION FREQUENCIES TO THE' 
CORRESPONDING DIGITAL PROTOTYPE FREQUENCIES ' 
NOTE: YOU WILL NEED PROTOTYPE FREQUENCY • 
RESPONSE PLOTS LIKE THOSE GIVEN IN REFC2) FOR ' 
BUTTERWORTH FILTERS C FIG 10.55, PGS 6918692) ' 

AND CHEBYSHEV FILTERS CFIG 10.54, PGS 696-696) ' 



READC 5 , X ) TYPE 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
CALL EXCMS C 'CLRSCRN') 

WR I T EC 6 , X ) 

WRITEC6,X) 1 WHAT IS THE ORDER OF THE PROTOTYPE YOU WISH TO USE?' 
WRITEC6,x) 

WRITEC 6 , X ) 

WRITEC6,X) 

WRITEC 6 , X ) 

WRITEC 6 , X) 

WRITEC 6 , X) 

WRITEC 6 , x) 

WRITEC 6 , X) 

WRITEC 6 , X) 

WRITEC 6 , x) 

WRITEC 6 , x) 

R EADC 5 , X ) ORDER 

IFC C TYPE. GT. 4) .AND. CTYPE.LT. 8) )THEN 
IFC ORDER . EQ . 1 )THEN 

WRITEC 6 , X) 'ELLIPTIC FILTERS GIVEN FOR 2ND-9TH ORDER ONLY' 

WRIT EC 6 , x) 

GO TO 10 
ENDIF 
ENDIF 

CALL EXCMS C 'CLRSCRN' ) 

MR I TEC 6 , x) 'DO YOU WANT TO ENTER SPECIFICATION FREQUENCIES' 
WRITEC 6 , x) 'AS ANALOG OR DIGITAL FREQUENCIES?' 

WRI TEC 6 , X) ' 1-ANALOG ' 

WRI TEC 6 , X) • 2-DIGITAL ' 

WRITEC6,*) ' XX INTEGER 1 OR 2 xx • 

R EADC 5 , x ) AFORDF 
CALL EXCMS C 'CLRSCRN' ) 

IFCAFORDF.EQ.DTHEN 

WRI TEC 6 , X ) 'WHAT IS THE SAMPLING FREQUENCY IN KHZ?' 

WRI TEC 6 , x ) « xx DECIMAL XX* 

READC 5 , x ) SSFREQ 
SFREQ=DBLECSSFREQ) 

ENDIF 

CALL EXCMS C 'CLRSCRN') 

WRI TEC 6, X) 

WRITEC 6 , x ) 

WRITEC 6 , x ) • 

WRI TEC 6 , X ) • 

WRI TEC 6 , x ) • 

WRITEC 6 , x ) • 

WRITEC6, X) 

WRITEC6 , X) • 

READC 5 , x) PASS 

FORMULAS FOR HIGHPASS AND LOWPASS 
IFC CPASS. EQ. 1) .OR. CPASS. EQ.2))THEN 
CALL EXCMS C 'CLRSCRN') 

IFCAFORDF.EQ.DTHEN 

WRI TEC 6 , X ) 'WHAT IS THE CRITICAL FREQUENCY IN KHZ? ' 

MR I T E C 6 , X ) ' - NEG 5DB PT FOR BUTTERWORTH ' 

WRITEC 6 , x) ' - RIPPLE EDGE FOR ELLIPTIC AND CHEBYSHEV' 

WRITE C 6 , X) ' XX DECIMAL XX ' 

READC 5 , x ) S AFC 
AFC r DBL EC SAFC) 

DFC= C AFC/S FREQ) XPIXC 2. 0D+00) 

ENDIF 

IFCAFORDF. EQ.2UHEN 

WRI T EC 6 , X) 'WHAT IS THE CRITICAL DIGITAL FREQUENCY?' 



'WHAT TYPE OF PASSBAND DUE YOU WANT?' 
' 1 -L OWPASS ' 

' 2-HIGHPASV 

' 5-BANDPASS' 

' 4-BANDSTOP' 

xx INTEGER 1-4 xx » 



NEW00750 
NEW00740 
NEW00750 
NEW0076 0 
NEW0077 0 
NEW0078 0 
NEW00790 
NEW00800 
NEW005 1 0 
NEW0082 0 
NEW0085 0 
NEW00840 
NEW008 5 0 
NEW00860 
NEW0087 0 
NEW00880 
NEW00890 
NEW00900 
NEW00910 
NEW0092 0 
NEW00950 
NEW0094 0 
NEW00950 
NEW00960 
NEW00970 
NEW00980 
NEW00990 
NEMO 10 00 
NEW01010 
N EM 0 1 0 2 0 
NEW01050 
NEW01040 
NEW01050 
NEW01060 
NEMO 107 0 
NEMO 1 080 
NEM01090 
NEW01100 
NEW01110 
NEW01120 
NEMO 1150 
NEW01140 
NEMO 1150 
NEMO 1 16 0 
NEMO 117 0 
NEW01180 
NEMO 119 0 
NEMO 120 0 
NEW012 1 0 
NEMO 1 220 
NEMO 1250 
NEMO 1 24 0 
NEMO 1 250 
NEMO 126 0 
NEMO 1 27 0 
NEMO 1280 
NEMO 1290 
NEMO 1 50 0 
NEMO 1 51 0 
NEW 01 520 
NEW01550 
NEW01540 
NEMO 1550 
NEMO 1560 
NEMO 1570 
NEMO 1 380 
NEMO 1 59 0 
NEMO 1400 
NEW0141 0 
NEW01420 
NEW01450 
NEW01440 
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HRITEC6,*) * - NEG 3DB FOR BUTTERWORTH » 

WRITEC6,X) » - RIPPLE EDGE FOR ELLIPTIC AND CHEBYSHEV* 

WRITEC6,X) 1 xx DECIMAL XX ' 

READ(5,X) SDFC 
DFC = DBL EC SDFC) 

ENDIF 

DFCD2=DFC/2.0D+00 

CALCULATE THE APPROPRIATE VALUES OF THE DESIGN CONSTANT ALPHA 
FROM TABLE 10.7 REFERENCE (2) AND REFERENCE (3) 

IFC PASS . EQ . 1 )THEN 

ALPHA=DSINCCPID4)-(DFCD2))/DSINCCPID4)+(DFCD2)) 

ENDIF 

IFCPASS.EQ.2)THEN 

AL PHA = DCOSC (PID4)-(DFCD2) )/DC0S( (PID4)+(DFCD2)) 

ENDIF 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
OPTION TO SHOW CALCULATED VALUE OF ALPHA 
CALL EXCMS C*CLRSCRN*) 

WRITEC6,X) * DO YOU WISH TO SEE THE CALCULATED VALUE OF ALPHA?* 
WRI TEC 6 , X) • XX 1 -YES/ 2-NO xx» 

READC 5 , X) LOOK 
IFC LOOK . EQ . 1 )THEN 
WRI TEC 6 , x ) 

WRI TEC 6 , 903 ) ALPHA 
WRITEC 6 , x) 

WRITEC6/X) *xx ENTER 1 TO CONTINUE XX* 

READ C 5 , X) LOOK 
ENDIF 
ENDIF 

xxxxxxxxxxxxx* rxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
BANDPASS AND BANDSTOP FORMULAS 
I FC C PASS . EQ . 3) . OR . C PASS . EQ . A ) ) THEN 
CALL EXCMS C ’CLRSCRN*) 

IFCAFORDF. EQ . 1 )THEN 

WRI TEC 6 , X ) * WHAT IS THE LOWER CRITICAL FREQUENCY IN KHZ?* 
WRITEC 6 / x) » - NEG 3DB FOR BUTTERWORTH* 

WRITEC 6 / X) * - RIPPLE EDGE FOR ELLIPTIC AND CHEBYSHEV » 

WRITEC6/X) » xx DECIMAL xx» 

READC 5 , x ) SALC 
ALC=DBL EC SALC ) 

CALL EXCMS C * CLRSCRN * ) 

WRITEC6/X) * WHAT IS THE UPPER CRITICAL FREQUENCY IN KHZ?* 

WRI T EC 6 / x ) * - NEG 3DB FOR BUTTERWORTH* 

WRITEC 6 , x ) » - RIPPLE EDGE FOR ELLIPTIC AND CHEBYSHEV » 

WRITEC6/X) * xx DECIMAL xx» 

READC 5# x ) SAUC 
AUC= DBL EC SAUC) 

DLC=CALC/SFREQ)XC2. 0D+00)XPI 
DUC r (AUC/SFREQ)x(2 . 0 D+ 00 )XPI 
ENDIF 

IFCAFORDF. EQ.2)THEN 

WRITEC 6 , x ) * WHAT IS THE LOWER CRITICAL DIGITAL FREQUENCY?* 
WRITEC 6 / x ) ' - NEG 3DB FOR BUTTERWORTH * 

WRITEC 6 , x ) » - RIPPLE EDGE FOR CHEBYSHEV AND ELLIPTIC* 

WRI TEC 6 , x) » xx DECIMAL xx • 

READC5/X) SDLC 
DLC = DBLEC SDLC) 

CALL EXCMS C'CLRSCRN*) 

WRITEC 6 , x ) » WHAT IS THE UPPER CRITICAL DIGITAL FREQUENCY?* 
WRITEC 6 , x ) » - NEG 3DB FOR BUTTERWORTH * 

WRITEC6/X) * - RIPPLE EDGE FOR CHEBYSHEV AND ELLIPTIC* 

WRITEC6 , x) » xx DECIMAL xx * 

READC 5 , x ) SDUC 
DUC = DBL EC SDUC ) 

ENDIF 

DLCD2=DLC/2.0D+00 

DUCD2=DUC/2.0D+00 

CALCULATE APPROPRIATE VALUES OF DESIGN CONSTANTS ALPHA AND K 
FROM TABLE 10.7 REFERENCE C2) AND REFERENCE (3) 
ALPHA=DC0SCCDUCD2)+CDLCD2))/DC0SCCDUCD2)-CDLCD2) ) 

IFC PASS . EQ . 4 )THEN 

RK=DTANCPID4)XDTAN( ( DUCD2 ) - C DLCD2 ) ) 



NEW01A50 
NEW01460 
NEW01A70 
NEW01480 
NEW01A90 
NEW01500 
NEW01510 
NEW01520 
NEW01530 
NEW015A0 
NEW01550 
NEW01560 
NEW01570 
NEW01580 
NEW01590 
NEW01600 
NEW01610 
NEW01620 
NEW01630 
NEW016A0 
NEW01650 
NEW01660 
NEW0167 0 
NEW01680 
NEW01690 
NEW01700 
NEMO 17 1 0 
NEW01720 
NEW01730 
NEW017A0 
NEW01750 
NEW01760 
NEW01770 
NEW01780 
NEW017 90 
NEW01800 
NEW01810 
NEW01820 
NEW01830 
NEW018A0 
NEW01850 
NEW01860 
NEW0187 0 
NEW01880 
NEW01890 
NEW01900 
NEW01910 
NEW01920 
NENO 1930 
NEW019A0 
NEW01950 
NEW01960 
ME WO 1970 
NEW01980 
NEW01990 
NEW02000 
NEW0201 0 
NEW02020 
NEW02030 
NEW020A0 
NEW02050 
NEW0206 0 
NEW02070 
NEW02080 
NEW02090 
NEW02100 
NEW02110 
NEW021 20 
NEW02130 
NEW021A0 
NEW02150 
NEW0216 0 



ii 



37 



oooo 



II 



FILE* DFCADD FORTRAN A 



C 



16 



B=( C2 . 0D+00)*ALPHA)/(RK+1 .0D+00) 
C=CC1.0D+00)-RK)/CC1.0D+00)+RK) 

ENDIF 

IFCPASS. EQ.3)THEN 

RK=DTAN(PIDA)/DTAN( ( DUCD2 )-C DLCD2 ) ) 

B = CC2 .0D+00)xALPHAxRK)/(RK+C1.0D+00)) 

C=CRK-C1.0D+00) )/CRK+(l . 0D+00) ) 

ENDIF 

OPTION TO SHOW CALCULATED VALUE OF ALPHA AND K 
CALL EXCMS C'CLRSCRN') 

WRITEC6,*) 'DO YOU WISH TO SEE THE CALCULATED VALUE OF ALPHA * 
WRITE<6,*) ' AND K? ' 

WRITEC 6 , x ) • xx 1-YES/2-N0 xx* 

READC5,x) LOOK 
IFCLOOK. EQ.DTHEN 
WRITEC6,*) 

WRI TEC 6 ,903) ALPHA 
WRI TEC 6 , 9 0 A ) RK 
WRITEC6 , X) 

WRIT EC 6 , x) 'xx ENTER 1 TO CONTINUE xx* 

READC 5 , X ) LOOK 
ENDIF 
ENDIF 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
OPTION TO ASSIST IN DETERMINING THE PROTOTYPE ORDER BY CHANGING 
SPECIFICATION FREQUENCIES TO APPROPRIATE PROTOTYPE FREQUENCIES 
FROM EQN 10.27A REFERENCEC 1 ) 

I FC ORDER . EQ . 0 )THEN 
CALL EXCMS C'CLRSCRN') 

IFCAFORDF. EQ.DTHEN 

WRITEC6 , X) 'WHAT IS THE ANALOG SPECIFICATION FREQUENCY (KHZ)' 
WRI TEC 6 , X ) ' THAT YOU WANT TO CONVERT TO A DIGITAL' 

WRI TEC 6 , X ) 'LOWPASS PROTOTYPE FREQUENCY?' 

WRI T EC 6 , x ) • xx DECIMAL Xx» 

READC 5 , X) F 

DF=C2.0)x$plx(F/S$FREQ) 

ENDIF 

IFCAFORDF. EQ.2)THEN 

WRI TEC 6 , X) ' WHAT IS THE DIGITAL SPECIFICATION FREQUENCY' 

WRI TEC 6 > X ) ' THAT YOU WANT TO CONVERT TO A DIGITAL' 

WRI TEC 6 , X) 'LOWPASS PROTOTYPE FREQUENCY?' 

WRI T EC 6 , x ) • xx DECIMAL xx» 

READC 5 , X ) DF 
ENDIF 
X=COSCDF) 

Y=SI N C DF ) 

Z = CMPL XC X/ Y) 

IFCPASS. EQ.DTHEN 
AS = REAL C ALPHA ) 

NEHZ=CZ-AS)/CC1 . 0)-CASXZ)) 

ENDIF 

IFCPASS. EQ.2)THEN 
AS = REAL C ALPHA ) 

NEWZ=C-1 .0)X(CZ-AS)/C Cl .0)-CASXZ))) 

ENDIF 

IFCPASS. EQ.3)THEN 
BS = REAL C B ) 

CS = REAL C C) 

NEWZ=C-1. 0)XCCCZXX2)-CBSXZ)+CS)/C Cl .0)-CBSXZ) + CCSXCZxx2)))) 
ENDIF 

IFCPASS. EQ. A)THEN 
BS = REAL C B) 

CS = REA L C C ) 

NEWZ=CCZXX2)-CBSXZ)+CS)/(C1 .0)-CBSxZ)+CCSXCZXX2)D 
ENDIF 

PTHETA = ABSC ATANC AIMAGC NEWZ)/REAL C NEWZ) ) ) 

IFCPASS. EQ.DTHEN 
IFC F . GT . SAFOTHEN 
PTHETA=SPI-PTHETA 
ENDIF 
ENDIF 

IFCPASS. EQ.2)THEN 



NEW0217 0 
NEW02180 
NEW02190 
NEW02200 
NEW0221 0 
NEW02220 
NEW02230 
NEW022A0 
NEW02250 
NEW02260 
NEW02270 
NEW02280 
NEW0229 0 
NEW02300 
NEW02310 
NEW02320 
NEW02330 
NEW023A0 
NEW02350 
NEW 02 36 0 
NEW0237 0 
NEW0238 0 
NEW02390 
NEW02A00 
NEW02A1 0 
NEW02A20 
NEW02A30 
NEW02AA0 
NEW02A50 
NEW02A6 0 
NEW02A70 
NEN02A80 
NEW02A 90 
NEW02500 
N EM 0 2 5 1 0 
NEW02520 
NEW02530 
NEN025A0 
NEW02550 
NEW 0 256 0 
NEW 02 570 
NEW 02 580 
NEW02590 
NEW 026 00 
NEW 026 10 
NEW 026 20 
NEW 026 30 
NEW026A0 
NEW 026 50 
NEW0266 0 
NEW02670 
NEW 026 80 
N EH 026 9 0 
NEW02700 
N EW 02710 
NEW 027 20 
NEW02730 
NEW027 A 0 
NEW02750 
NEN0276 0 
NEW0277 0 
NEW0278 0 
NEW 027 90 
NEW028 0 0 
NEN02810 
NEW 02 820 
NEW 028 30 
NEW028A0 
NEW02850 
NEW02860 
NEW0287 0 
NEW02880 
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C 



IF( F . LT . SAFOTHEN 
PTHETA=SPI-PTHETA 
ENDIF 
ENDIF 

IF(PASS.EQ.3)THEN 

IF( (F.LT.SALC) . OR . C F . GT . SAUC) ) THEN 
PTHETA=SPI-P THETA 
ENDIF 
ENDIF 

IF( PASS . EQ . 4)THEN 

IF( (F.GT.SALC) . AND.CF.LT. SAUC))THEN 
PTH ETA=SPI-P THETA 
ENDIF 
ENDIF 

WRITEC6,905) PTHETA 
WRIT EC 6 , X ) 

WRIT EC 6 f X ) 'WANT TO CONVERT ANOTHER FREQUENCY?' 

WRITE C 6 , X ) ' XX 1-YES/2-N0 XX* 

READC 5 / X ) LOOK 
IFC LOOK . EQ . 1 )GO TO 16 
IFC LOOK . EQ . 2)THEN 
WR IT EC 6 , X ) 

WRIT EC 6 , X ) 1 WHAT ORDER PROTOTYPE DO YOU WANT TO USE?' 

WRITE C 6 , X ) ■ XX INTEGER 1-5 XX * 

READC5/X) ORDER 
ENDIF 
ENDIF 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
LOAD APROT MATRIX WITH APPROPRIATE COEFFICIENTS 
APROT C TYPE- 1 ) ARE BUTTERWORTH FILTERS 
IFC TYPE . EQ . 1 )THEN 
IFC ORDER . EQ . 1 )THEN 
APROT C 1 , 0 ) = 1 .0D+00 
APROT C 2 , 0 ) = 1 . 0 D+00 
APROT C 2 / 1 ) = 1 . 0 D+ 0 0 
ELSEIFC0RDER.EQ.2)THEN 
APROTC1,0)=1 .0D+00 
APROT C 2/ 0 ) = 1 . 0 D+0 0 
APR0TC2/1 ) = DSQRTC2. OD+OO) 

APR0TC2/2)=1 .OD+OO 
ELSEIFCORDER.EQ.3)THEN 
APR0TC1,0)=1. OD+OO 
APROT C 2 , 0 ) = 1 .OD+OO 
APR0TC2, 1 )=2. OD+OO 
APR0TC2/2)=2. OD+OO 
APROT C 2/ 3 ) = 1 .OD+OO 
ELSEIFCORDER.EQ.4)THEN 
APROT C 1 , 0 ) = 1 .OD+OO 
APR0TC2, 0) = 1 . OD+OO 
APROTC2,1)=2.613126D+00 
APR0TC2,2)=3. 414214 D+00 
APR0TC2,3)=2. 613126 D+00 
APROT C 2, 4 ) = 1 .OD+OO 
ELSEIFC0RDER.EQ.5)THEN 
APROT C 1 / 0 ) = 1 .OD+OO 
APROT C 2/ 0 ) = 1 .OD+OO 
APR0TC2,1 )=3. 236068 D+00 
APR0TC2,2)-5. 236068 D+00 
APROT C 2, 3) =5. 236 06 8 D+00 
APR0TC2,4)=3. 236068 D+00 
APROT C 2, 5 ) = 1 .OD+OO 
ENDIF 
ENDIF 

APROT C TYPE = 2 ) ARE CHEBYSHEV FILTERS WITH . 5DB RIPPLE 
I FC TYPE . EQ . 2 )THEN 
IFC ORDER . EQ . 1 )THEN 
APROTCl,0)-2. 862775 D+00 
APROTC2,0)=2. 862775 D+00 
APROT C 2/ 1 ) = 1 .OD+OO 
EL SEIFC ORDER . EQ . 2 )THEH 
APR0TC1,0)=1 .431388D+00 
APROT C 2 / 0 ) = 1 . 51 6 20 3D+0 0 



NEW02890 
NEW029 0 0 
NEW029 1 0 
NEW02920 
NEW02930 
NEW02940 
NEW02950 
NEW0296 0 
NEW0297 0 
NEW0298 0 
NEW0299 0 
NEW03000 
NEW03010 
NEW03020 
NEW03030 
NEW03 04 0 
NEW03050 
NEW03060 
NEW03070 
NEW0308 0 
NEW03090 
NEW03100 
NEW03110 
NEW03120 
NEW031 30 
NEW03140 
NEW031 5 0 
NEW03160 
NEMO 31 7 0 
NEW03180 
NEW03190 
NEN03200 
NEN03210 
NEW03220 
NEW03230 
NEW0324 0 
NEW03250 
NEW 0326 0 
NEW03270 
NEW0328 0 
NEW03290 
NEMO 330 0 
NEW03310 
NEMO 3320 
NEMO 3330 
NEW0334 0 
NEMO 3350 
NEN03360 
NEW03370 
NEW0338 0 
NEW03390 
NEMO 34 0 0 
NEW03410 
NEN03420 
N EM 034 3 0 
NEM0344 0 
NEW03450 
NEH0346 0 
NEW03470 
NEMO 3 48 0 
NEM03490 
NEMO 350 0 
NEW03510 
NEMO 3520 
NEN03530 
NEW03540 
NEW03550 
NEW0356 0 
NEMO 357 0 
NEM0358 0 
NEW03590 
NEW036 0 0 
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APRDTC2,l)=1.425625D+00 
APROT (2, 2) = 1 . 0D+00 
ELSEIFCORDER.EQ.3)THEN 
APROTC1,0)=0. 715694 D+ 00 
APR OTC2,0)=0. 715694 P+00 
APR0TC2,1)=1 . 534895D+00 
APROT ( 2, 2)=1 . 252913 D+ 00 
APR0TC2,3)=1 .0D+00 
ELSEIFC0RDER.EQ.4)THEN 
APR0TC1,0)=0. 357847 D+ 00 
APR0TC2,0)=0. 379051 D+ 00 
APROTC 2, 1 ) = 1 . 025455D+00 
APR0TC2,2)=1 .716866 D+ 00 
APR0TC2,3)=1 .197386D+00 
APR0TC2>4)=1 .0D+00 
ELSEIFC ORDER . EQ . 5 )THEN 
APROT ( 1 / 0 )=0 . 178923D+00 
APR0TC2, 0)=0.178923D+00 
APROTC 2, l)=0.752518D+00 
APR0T(2/2)=1 .309575 D+ 00 
APR0T(2,3)=1 .937368D+00 
APROTC 2, 4 ) = 1 . 172491D+00 
APR0T(2/5)=1 .0D+00 
ENDIF 
EfJDIF 

C APROT C T YPE=3 ) ARE CHEBYSHEV FILTERS WITH 1DB RIPPLE 

I F ( TYPE . EQ . 3 )THEN 
I F( ORDER . EQ . 1 )THEN 
APROTC 1 ^ 0 ) = 1 . 1 96523D+00 
APROTC 2, 0) = 1 . 1 96 523 D+ 00 
APR0TC2/1)=1. OD+OO 
ELSEIFC ORDER , EQ . 2)THEN 
APR0TCl,0)=0.982613D+00 
APROTC2,0)=1. 102510 D+ 00 
APR0TC2,1)=1 .097734 D+ 00 
APROT ( 2, 2)=1 .OD+OO 
ELSEIFC ORDER . EQ . 3)THEN 
APROTC 1, 0)=0. 491307 D+ 00 
APR 0TC2,0)=0. 491307 D+ 00 
APROTC2>1)=1.238409D+00 
APR0TC2,2)=0. 988341 D+ 00 
APROT C 2, 3 ) = 1 . 0D+00 
ELSEIFC ORDER . EQ . 4 )THEN 
APROTC 1 , 0 ) = 0 .245653D+00 
APR0TC2, 0)=0.275628D+00 
APROTC 2, l)=0.7426l9D+00 
APR0TC2,2)=1. 453925 D+ 00 
APROTC2,3)=0.952811D+00 
APR0TC2/4)=1 . OD+OO 
ELSEIFC ORDER . EQ.5)THEN 
APROTC1,0)=0. 122827D+00 
APROTC2,0)=0.122827D+00 
APROTC2,1)=0. 580534 D+ 00 
APROT C2,2)=0.974396D+00 
APR0TC2,3)=1 . 688816 D+ 00 
APROTC2,4)=0. 936820 D+ 00 
APR0TC2, 5)=1 . OD+OO 
ENDIF 
ENDIF 

C APROT C TYPE = 4 ) ARE CHEBYSHEV FILTERS WITH 2DB RIPPLE 

I F C TYPE . EQ . 4 )THEN 
IFC ORDER . EQ . 1 )THEN 
APROT C 1 , 0) = 1 . 307560 D+ 00 
APR 0TC2;0)=1. 307560 D+ 00 
APROT C2/l)-l. OD+OO 
ELSEIFC0RDER.EQ.2)THEN 
APROTC 1 > 0 ) = 0 . 505803 D+ 00 
APR0TC2, 0)=0 .636768D+00 
APROTC 2, 1 ) = 0 .803816D+00 
APROT C 2, 2) = 1 . OD+OO 
ELSEIFC ORDER . EQ . 3)THEN 
APROT C 1 , 0 ) = 0 .326890D+00 
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APROT(2,0)=0.326890D+00 
APR0T(2,1)-1 .022190D+00 
APROT(2,2)=0.737822D+00 
APR0T(2,3)-1 .0D+00 
ELSEIF(0RDER.EQ.4)THEN 
APR0T(l,0)=0.163445D+00 
APROT ( 2, 0)=0. 205765 D+ 00 
APR0T(2,1)=0. 516798 D+ 00 
APR0T(2,2)=1 .256482D+00 
APROT(2,3)=0.716215D+00 
APROT ( 2/ A ) =1 . 0 D+ 0 0 
ELSEIF(0RDER.EQ.5)THEN 
APROT ( 1 , 0 ) = 0 . 081723 D+ 00 
APROT(2,0)=0. 081723 D+ 00 ' 

APROT(2,1)=0. 459349 D+ 00 
APR0T(2,2)=0.693477D+00 
APR0T(2,3)=1 .499543D+00 
APR0T(2,4)=0.706461D+Q0 
APR0T(2,5)=1 .CD+00 
ENDIF 
ENDIF 

APROT ( TYPE = 5 ) ARE ELLIPTIC FILTERS 

ELLIPTIC COEFFICIENTS ARE CALCULATED USING A PROGRAM DEVELOPED 

BY PROFESSOR D. E. KIRK, NAVAL POSTGRADUATE SCHOOL, MONTEREY, 

CALIFORNIA 93943 

A02=0.0D+00 

B 02 = 0 . 0D+00 

B 1 2 = 0 . 0 D+ 0 0 

S0=0.0D+00 

IFC TYPE. EQ . 5) THEN 

CALL EXCMS C »CLRSCRN» ) 

WRITEC 6 , X) * WHAT PASSBAND RIPPLE DO YOU WANT? (IN DB)' 

WRITEC 6 , X) *xxx MUST BE EITHER 0.5, 1.0, 2.0, OR 3.0 ***• 
READC 5, X) ASUBP 
WRITEC 6, X) 

WRIT E( 6 , X ) 

WRIT EC 6 , x ) * WHAT STOPBAND ATTENUATION DO YOU WANT? (IN DB)» 
NRITEC 6 , X ) *xxx MUST BE “20,-30,-40,-50,-60, OR -70 XXX* 
WRITEC 6 , X ) »xx* DO NOT INCLUDE A DECIMAL XXX* 

READC 5, X ) NATTN 
NATTN=-1 XNATTN 
ASUBA=REAL ( NATTN) 

IFCASUBP .EQ.0.5)THEN 
IFC NATTN . EQ . 20)THEN 

IFCORDER. EQ.2)SR=2. 76261 
IFCORDER. EQ.3)SR=1 .42189 
IFCORDER. EQ.4)SR=1 .13188 
IFCORDER. EQ.5)SR=1 .04465 
IFCORDER. EQ.6)SR=1 .01553 
ENDIF 

IFC NATTN. EQ . 30)THEN 

I F( ORDER. EQ. 2)SR=4 .80880 
IF(0RDER.EQ.3)SR=1 .92322 
IFCORDER. EQ.4)SR=1 .32446 
IF(0RDER.EQ.5)SR=1 . 12912 
IFCORDER. EQ.6)SR=1 . 05394 
ENDIF 

IFC NATTN . EQ . 40 )THEN 

IFCORDER. EQ. 2 )SR=8 .848925 
IFCORDER. EQ.3)SR=2. 71147 
IFCORDER. EQ.4)SR=1. 62842 
IFCORDER. EQ.5)SR=1 .27264 
IFCORDER. EQ.6)SR=1 .12697 
ENDIF 

IFC NATTN. EQ. 50 )THEN 

IFCORDER. EQ. 2) SR= 15 . 06069 
IFCORDER. EQ.3)SR=3. 90430 
IFCORDER. EQ.4)SR=2. 06924 
IFCORDER. EQ.5)SR=1. 48469 
IFCORDER. EQ.6)SR=1 .24101 
ENDIF 

IFC NATTN . EQ .6 0)THEN 
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END1F 

1F(NATTN. EQ.40)THEN 

IFCORDER. EQ.2)SR=5. 76107 
IFCORDER. EQ.3)SR=2. 13923 
IF( ORDER .EQ.4)SR=1. 40842 
IFCORDER. EQ.5)SR=1 .16811 
IFCORDER. EQ.6)SR=1 .07316 
ENDIF 

1 F( NATTN . EQ . 50 ) THEN 

IFCORDER. EQ.2)SR=10. 19181 
IFCORDER. EQ.3)SR=3. 04137 
IFCORDER. EQ.4)SR=1 .75285 
IFCORDER. EQ.5)SR=1 .33243 
1FC0RDER.EQ.6)SR=1. 15868 
ENDIF 

I F C NATTN .EQ.60)THEN 

IFCORDER. EQ.2)SR=18 .09398 
IFCORDER. EQ.3)SR=4. 39729 
IFCORDER. EQ.4)SR=2. 24440 
IFCORDER. EQ.5)SR=1 .56860 
IFCORDER. EQ.6)SR=1 .28693 
ENDIF 

I F( NATTN. EQ . 70 )THEN 

IFCORDER. EQ.2)SR=32. 15951 
1FCORDER. EQ.3)SR=6 .40894 
IFCORDER. EQ.4)SR=2. 92363 
IFCORDER. EQ.5)SR=1 .88906 
IFCORDER. EQ.6)SR=1 .46330 
ENDIF 
ENDIF 

IFCASUBP.EQ.3.0)THEN 
I FC NATTN . EQ . 20 )THEN 

IFCORDER. EQ.2)SR=1 .73915 
1F(0RDER.EQ.3)SR=1. 15516 
IFCORDER. EQ.4)SR=1. 03853 
IF( ORDER . EQ . 5 ) SR=1 .00996 
IFCORDER. EQ.6)SR=1 .00260 
ENDIF 

I F( NATTN . EQ . 30 )THEN 

IFCORDER. EQ.2)SR=2. 90352 
IFCORDER. EQ.3)SR=1. 45814 
IFCORDER. EQ.4)SR=1. 14542 
IFCORDER. EQ.5)SR=1 .05020 
IFCORDER. EQ.6)SR»1 .01783 
ENDIF 

1F( NATTN . EQ . 40) THEN 

IFCORDER. EQ.2)SR=5. 05584 
IFCORDER. EQ.3)SR=1 .98022 
IFCORDER. EQ.4)SR=1 .34663 
IFCORDER. EQ.5)SR=1 .13934 
IFCORDER. EQ.6)SR=1 .05892 
ENDIF 

1 F( NATTN . EQ . 50 ) THEN 

IFCORDER. EQ.2)SR=8. 93009 
IFCORDER. EQ.3)SR=2. 79862 
IFCORDER. EQ.4)SR=1 .66148 
IFCORDER. EQ.5)SR=1 .28850 
IFCORDER. EQ.6)SR=1 .13534 
ENDIF 

I F( NATTN . EQ . 60)THEN 

IFCORDER. EQ.2)SR=15. 84610 
IFCORDER. EQ.3)SR=4. 03471 
IFCORDER. EQ.4)SR=2. 11596 
IFCORDER. EQ.5)SR=1 .50711 
IFCORDER. EQ.6)SR=1 .25325 
ENDIF 

1 FC NATTN . EQ . 7 0 ) THEN 

IFCORDER. EQ.2)SR=28. 15950 
IFCORDER. EQ.3)SR=5. 87251 
IFCORDER. EQ.4)SR=2. 74749 
IFCORDER. EQ.5)SR=1 .80680 
IFCORDER. EQ.6)SR=1 .41799 
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17 



19 



ENDIF 

ENDIF 

SPI =3 .14159 
RK=1 .0/SR 
N=ORDER 
EN = N 

KPR=SQRT(1.0-RK**2) 

Q0=0.5X(1 .0-SQRT(KPR))/(l . Q + SQRT ( KPR ) ) 

Q = Q0+2 . 0*Q0x*5+15. 0XQOXX9+150 . OXQOXX13 
D=(10.0XX(0 .1XASUBA)-1.0)/(10.0XX(0. 1XASUBP)-1 .0) 

LAM= ( 1 . 0/(2.0XEN))XALOG((10.0XX(0.05*ASUBP)+1 . 0 )/( 10 . Oxx( 0 . 05 

C XASUBP ) “1 . 0 ) ) 

I STOP = 0 
DENSUM=0 .0 
ENSUM=SINH( LAM) 

M = 1 
EM=M 

ENUM=((-1 .0)X*M)X(QXX(MX(M+1)))XSINH((2.0XEM+1 .0)*LAM) 
DEN= 2 . 0 X( (-1 .0)x*M)x(Qxx(Mxx2))xCOSH(2.0xEMxLAM) 
ENSUM=ENSUM+ENUM 
DENSUM=DENSUM+DEN 
I F( M . GE . 2)THEN 

RN=ABS( ENUM/ENSUM) 

RD=ABS(DEN/(1.0+DENSUM)) 

I F( ( RN . LE . 1 .OE-8 ).AND.(RD.LE.l. OE-8 ) )THEN 
ISTOP=l 
ENDIF 
ENDIF 
M=M+1 

IF(ISTOP.EQ.O)THEN 
GO TO 17 
ENDIF 

SIGMAO=ABS((2.0X(QX*0.25)XENSUM)/C1 .0+DENSUM)) 

W=SQRT(C1 .0+RKX(SIGMAOXX2))X(1 .0+CSIGMA0XX2/RK))) 

IN=M0D(N,2) 

IF( IN . EQ . 0 )THEN 
IR=N/2 
ELSE 

IR=(N-l)/2 

ENDIF 

DO 18 1=1, IR 

I F( IN . EQ . 0 )THEN 
EMU=I-0 . 5 
ELSE 
EMU = I 
ENDIF 
I STOP=0 
DENSUM=0 . 0 

ENSUM=SIN( PixEMU/EN ) 

M=1 
EM = M 

ENUM= ( ( -I . 0)XXM)X(QXX(MX(M+1)))XSIN((2.0XEM+1.0)XPIXEMU/EN) 

DEN=2. 0X((-1 . 0)XXM)X(QXX(MXX2))XCOS(2.0XEMXPIXEMU/EN) 

ENSUM=EIISUM+ENUM 

DENSUM=DENSUM+DEN 

IF( M.GE . 2)THEN 

RH=ABS( ENUM/ENSUM) • 

RD=ABS ( DEN/( 1 • 0+DENSUM) ) 

IF( CRN .LE.l .OE-8). AND. (RD.LE.l .OE-8) )THEN 
IST0P=1 
ENDIF 
ENDIF 
M r M+l 

IFC ISTOP.EQ.Q)THEN 
GO TO 19 
ENDIF 

OMEGA ( I)=(2.0X(Qxx0 . 25)xENSUM)/( 1 . O+DENSUM) 

V(I)=SQRT((1 .0-RKXOMEGA(I)XX2)X(1 . 0 - ( OMEGAC I ) XX2/RK) ) ) 

AOC I ) = 1 . 0/(0 MEG AC I )XX2) 

DN — 1 . 0+(SIGMAOXOMEGA(I))XX2 

B0(I)=((SIGMA0XVCI))XX2+C0MEGA(I)XW)XX2)/(DNXX2) 

B1(I)=(2.0XSIGMAOXVCI))/DN 



NEW06490 
NEH0650 0 
NEW06510 
NEW06520 
NEW06530 
NEW06540 
NEW06550 
NEW06560 
NEH0657 0 
NEW06580 
NEW06590 
NEW06600 
NEW06610 
NEH06620 
NEW06630 
NEW06640 
NEW06650 
NEW0666 0 
NEW0667 0 
NEW06680 
NEW06690 
NEH067 0 0 
NEN06710 
NEW06720 
NEWO6730 
NEW067A0 
NEW067 5 0 
NEN06760 
NEW0677 0 
NEW06780 
NEW06790 
NEH06800 
NEH06810 
NEH06820 
NEH06830 
NEW068A0 
N EN 06850 
NEW0686 0 
NEW0687 0 
NEW0688 0 
NEW06890 
NEW06900 
NEN06910 
NEW06920 
NEH06930 
NEW06940 
NEW06950 
NEH0696 0 
NEM06970 
NEH06980 
NEN06990 
NEW07000 
NEW07010 
NEN07020 
NEH07030 
NEH070AO 
NEH07050 
NEH07060 
NEMO 7 07 0 
NEM07080 
NEH07090 
NEH071 00 
NEW07110 
NEM07120 
NEH07130 
NEH071A0 
NEW07150 
NEW07160 
NEW07170 
NEW07180 
NEW07190 
NEW07200 



44 



FILEr DFCADD FORTRAN A 



18 



20 



21 



CONTINUE 
HO= 1 . 0 

IF( IN . EQ . 0 )THEN 
DO 20 1 = 1, IR 

HO=HO*BOCI)/AOCI) 

CONTINUE 

HO=HO*C10. 0**C-0 . 05*ASUBP)> 

ELSE 

DO 21 1=1, IR 

HO=HO*BOCI)/AOCI) 

CONTINUE 
HO=HO#SI GMAO 
ENDIF 

R = DBL E( SR) 

HO = DBL E( HO ) 

AO 1 = DB L E( AO ( 1 ).) 

B01 = DBL E( B0( 1 ) ) 

B1 1 =DBLE( B1 ( 1 ) ) 

IFCC0RDER.EQ.2) .OR. C ORDER . EQ . 3 )) THEN 
APROTCl,0)=H0*A01 
APROTC 1 , 1 )=0 . OD+OO 
APROT ( 1 , 2 ) =H0 
APROT (1,3)=0. OD+OO 
I F( ORDER . EQ . 2)THEN 
APROT (2,0)=B01 
APROT (2,1 )=B1 1 
APR0T(2,2)=1 .OD+OO 
ENDIF 

I F ( ORDER . EQ . 3) THEN 
S0 = DBL EC SIGMAO) 

APROTC2,0)=S0*B01 
APR0TC2, l)=B01+( B1 1 #S0 ) 

APROT(2,2)=SO+B11 
APR0T(2,3)=1 .OD+OO 
ENDIF 
ENDIF 

IF( (ORDER. EQ. A) .OR. C ORDER . EQ . 5 )) THEN 
A02=DBLE(AO(2)) 

B02=DBLE(BO(2)) 

B12=DBLE(B1(2)) 

APROTCl,0)=H0*A01*A02 
APROT(1,1)=0. OD+OO 
APROTCl,2)=H0*CA01+A02) 

APROTC1,3)=0. OD+OO 
APROT ( 1 , A ) =H0 
APR0T(1,5)=0. OD+OO 
I F ( ORDER . EQ . A)THEN 
APROTC2,0)=B01*B02 
APROTC 2, 1 )=( Bll*B02)+CB01*B12) 

APROT C2,2)=B01+B02+CB11*B12) 

APR0T(2,3)=B11+B12 
APROT ( 2 , A) =1 . OD+OO 
ENDIF 

I F( ORDER . EQ . 5 )THEN 
S0 = DBL EC SIGMAO ) 

APROTC2,0)=S0*B01*B02 

APROTC 2, 1 ) = C SOKC C B1 1 * B02 ) + C B 0 1 * B 12 ) ) ) + C BO 1 *B 02 ) 
APROTC2,2)=CS0*CB01+B02+CBll*B12)))+CBll*B02)+CB01*B12) 
APROTC 2, 3)=CS0*CB11+B12))+B01+B02+CB11*B12) 
APR0TC2,A)=SO+B11+B12 
APR0TC2,5)=1. OD+OO 
ENDIF 
ENDIF 

IFC C ORDER . EQ . 6 ) . OR . C ORDER . EQ . 7 ) )THEN 
A02=DBLECAOC2)) 

A03 = DBLECAOC3)) i,* 

B02=DBLECBOC2)) 

B03=DBLECB0C3)) 

B12=DBLECB1C2)) 

B13=DBLECB1C3)) 

APROTC 1,0 )=HO*A01*A02*A03 
APR0TC1,1)=0. OD+OO 



NEN07210 
NEH07220 
NEW07230 
NEW072A0 
NEW07250 
NEW07260 
NEW07270 
NEN07280 
NEW07290 
NEH07300 
NEN07310 
NEH07320 
NEW07330 
NEW073A0 
NEW07350 
NEH07360 
NEW07370 
NEW07380 
NEW07390 
NEW07A00 
NEW07A10 
NEN07A20 
NEW07A30 
NEN07AA0 
NEN07A50 
NEH07A60 
NEH07A70 
NEH07A80 
NEW07A90 
NEW07500 
NEH0751 0 
NEH07520 
NEW07530 
NEN075A0 
NEH07550 
NEH07560 
NEN07570 
NEW07580 
NEH07590 
NEH07600 
NEW07610 
NEH07620 
NEH076 30 
NEMO 76 AO 
NEN07650 
NEH0766 0 
NEW07670 
NEMO 76 8 0 
NEH07690 
NEH077 00 
NEH07710 
NEH07720 
NEH07730 
NEN077A0 
NEH07750 
NEW07760 
NEN077 7 0 
NEH07780 
NEN07790 
N EN078 0 0 
NEW07810 
NEN07820 
NEW07830 
NEN078A0 
NEH07850 
NEN0786 0 
NEN0787 0 
NEH0788 0 
NEH0789 0 
NEN07900 
NEN07910 
NEW07920 
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C 



26 

25 



61 

60 



APR0T(1,2)=H0*( (A02*A03)+(A01*CA02+A03))) 

APR0TCl/3)=0. OD+OO 

APROT (1, A)=HOX( A01+A02+A03) 

APROT<1,5)=0 .OD+OO 
APROT ( 1 , 6 ) =H0 

I FC (ORDER. EQ.6) .OR. (ORDER. EQ.7)) THEN ' 

APR0TC2,0)=B01XB02XB03 

APROT (2, 1)=(B11#B02*B03)+(B01*((B02*B13)+(B12*B03))) 
APROTC2,2)=CB02XB03)+CB11X((B02XB13)+CB12XB03)))+ 

C CB01X((B02+B03)+CB12XB13))) 

APROTC2/3)=B01XCB12+B13)+CB11XCCB02+B03)+CB12XB13)))+ 

C C(B02*B13)+(B12XB03)) 

APROT (2, A)=B01+CB11XCB12+B13))+C C B02+B03)+C B1 2XB1 3) ) 
APR0T(2,5)=B11+B12+B13 
APR0TC2,6)=1. OD+OO 
ENDIF 

IFC0RDER.EQ.7)THEN 
SO=DBLE( SIGMAO ) 

APR0TC2,7)=1. OD+OO 

APROTC2/6)=CSOXAPROTC2,6))+APROTC2,5) 
APROTC2,5)=CSOXAPROTC2,5))+APROTC2,A) 
APROTC2/A)=CSOXAPROTC2/A))+APROTC2,3) 
APROTC2,3)=CSOXAPROTC2/3))+APROTC2,2) 
APR0TC2/2)=CS0XAPR0TC2,2))+APR0TC2, 1) 

APROT (2, 1)=CSOXAPROTC2, 1))+APR0T(2, 0) 
APROTC2,0)=SOXAPROTC2/ 0) 

ENDIF 

ENDIF 

ENDIF 

NORMALIZE THE ELLIPTIC PROTOTYPES 
IFC TYPE . EQ . 5)THEN 
DO 25 1=1,2 
DO 26 J=0, ORDER 

APROTCI, J)=APROTCI, J)/( C DSQRT C R ) ) XX J ) 

CONTINUE 

CONTINUE 

ENDIF 

APROT ( TYPE=6 ) ALLOWS THE INPUT OF ANY TYPE OF PROTOTYPE FILTER 
BY IMPUTING THE COEFFICIENTS OF S 
I FCTYPE. EQ .6)THEN 
HR IT EC 6 , X) 

WRITEC6,*) 'THIS OPTION ALLOWS YOU TO USE ANY ANALOG ' 

WRITEC 6 , X) ' PROTOTYPE OF ORDER 1-10 1 

WRITEC 6 , x ) »xx COEFFICIENTS MUST BE IN DOUBLE PRECISION ' 

WRIT EC 6 , x ) * FORMAT C 0 . 123A56789D+00 ) XX ■ 

WRITEC6 , X) 

HRITEC 6 , X ) 'INPUT THE NUMERATOR COEFFICIENTS FIRST' 

DO 50 1=0/ ORDER 

WRITEC 6 , X ) * WHAT IS THE COEFFICIENT OF SXX',1,'?' 

READC5/X) APROT C 1 , 1 ) 

CONTINUE , 

WRITEC6 / X) 

CALL EXCMS C'CLRSCRN') 

WRITEC 6 , x ) 1 NOW INPUT THE DENOMINATOR COEFFICIENTS' 

WRITEC 6 / X ) 

DO 51 J = 0 / ORDER 

WRITEC 6 , X) 'WHAT IS THE COEFFICIENT OF SXX',J,'?» 

READC5/X) APROT C 2, J ) 

CONTINUE 

ENDIF 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
DEFINE ELEMENTS OF THE BILINEAR TRANSFORMATION MATRIX 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
FOR ALL THE BILINEAR TRANSFORATION MATRICES, THE FIRST COLUMN 
HAS ALTERNATING +1'S AND -l'S. THE LAST COLUMN HAS ALL +1'S 
DO 60 1=1/ ORDER 
DO 61 J=0/I 

BILINRCI, J,0)=(-1 .0)XXJ 
BILINRCI, J,I)=1.0 
CONTINUE 
CONTINUE 



NEW07930 
NEW079AO 
NEW07950 
NEW07960 
NEW07970 
NEW07980 
NEW07990 
NEW08000 
NEW0801 0 
NEW08020 
NEW08030 
NEW080AO 
NEW08050 
NEW08060 
NEW08070 
NEW 08 080 
NEW08090 
NEW08100 
NEW0811 0 
NEW08120 
NEW081 30 
NEW081 AO 
NEW08150 
NEW08160 
NEW08170 
NEW08180 
NEW08190 
NEW08200 
NEW0821 0 
NEW08220 
NEW08230 
NEW082AO 
NEW08250 
NEW08260 
NEW0827 0 
NEW08280 
NEW08290 
NEW08300 
NEW0831 0 
NEW08320 
NEW08330 
NEW083A0 
NEW08350 
NEW0836 0 
NEW0837 0 
NEW 08 380 
NEW 08 390 
NEW08A00 
NEW0841 0 
NEW08A20 
N EW08A30 
NEW08AA0 
NEW08A50 
NEW08A6 0 
NEW08A7 0 
NEW08A80 
NEH08A90 
NEW08500 
NEW0851 0 
NEW08520 
NEW08530 
NEW085A0 
NEH08550 
NEW0856 0 
NEW0857 0 
NEW08580 
NEW08590 
NEW086 00 
NEW086 1 0 
NEW086 20 
NEW08630 
NEW086A0 
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C 

C 

C 



C 



77 

76 

C 



78 

75 

C 

C 



81 

80 

C 

C 

C 



102 

101 

100 

C 

C 



C 



106 

105 



111 



112 



C 

C 

C 

C 

C 

C 



FIRST ORDER MATRIX HAS ALREADY BEEN LOADED CONTINUE INTERATION 
UNTIL YOU HAVE THE APPROPRIATE MATRIX 
DO 75 K=2, ORDER 
KM1 =K-1 

LOAD ELEMENTS IN INTERIOR COLUMNS, ALL ROWS EXCEPT LAST ONE 
DO 76 L=0,KM1 
DO 77 M = 1 , KM1 
MM1=M-1 

BILINRCK/L, M) =BI L INRC KM1 , L , M) + BI L INR C KM1 , L ,MM1 ) 

CONTINUE 

CONTINUE 

NOW LOAD THE LAST ROW 
DO 78 N = 1 f KM1 

BILINRCK,K,N)=BILINRCK,0,N)xCjC-1.0D+00)XX(K+N)) 

CONTINUE 

CONTINUE 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
THE FOLLOWING CLEARS DPROT PRIOR TO LOADING 
DO 80 K = 1 , 2 
DO 81 L = 0 , 1 0 
DPROT C K , L ) = 0 . 0 
CONTINUE 
CONTINUE 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
DIGITIZE THE ANALOG PROTOTYPE USING THE BILINEAR TRANSFORMATION 
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
DO 100 N = 1 , 2 
DO 101 J = 0, ORDER 
SUMR= 0 . 0 

DO 102 1=0, ORDER 

TEMP=APROT(N, I )*BI L INR( ORDER, I , J) 

SUMR=SUMR+TEMP 
CONTINUE 
DPROT C N, J ) =SUMR 
CONTINUE 
CONTINUE 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
OPTION TO SHOW DIGITAL PROTOTYPE COEFFICIENTS 
CALL EXCMS ( * CL RSCRN * ) 

WRITEC6,*) * THE DIGITAL PROTOTYPE HAS BEEN COMPUTED. DO YOU 1 
WRITEC6,*) ’WISH TO CHECK THE PROTOTYPE COEFFICIENTS?’ 

WRITEC 6 , X ) ’ XX 1-YES/2-N0 XX* 

READC5,X) LOOK 
I F( LOOK . EQ . 1 ) THEN 

MAKE DENOMINATOR COEFFICIENT OF ZXXORDER “1.0 
DO 105 K = l,2 
DO 106 L = 0, ORDER 

DPROTCK,L )=DPR0TCK,L)/DPR0TC2,0RDER) 

CONTINUE 
CONTI NUE 

CALL EXCMS (’CLRSCRN*) 

WRITEC 6 , x ) ’THE DIGITAL PROTOTYPE NUMERATOR COEFFICIENTS ARE t • 

DO 111 N=ORDER, 0,-1 
WRI T EC 6 , 9 02 ) DPR0TC1,N),N 
CONTINUE 
WRITEC 6 , X) 

WRIT EC 6 , X ) * THE DIGITAL PROTOTYPE DENOMINATOR COEFFICIENTS ARE j ' 
DO 112 M=ORDER ,0,-1 
WRITEC 6 # 902) DPR0TC2,M),M 
CONTINUE 
WRITEC 6 , x) 

WRITEC6,x) * x x 1 TO CONTINUE xx» 

READC 5, x) LOOK 
ENDIF 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
NOW DO THE DIGITAL TO DIGITAL TRANSFORMATIONS 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
IFC CPASS. EQ. 1 ) . OR . C PASS . EQ . 2) )THEN 

DEFINE ELEMENTS OF THE LOWPASS AND HIGHPASS DIGITA L -DIGI TAL 
TRANSFORMATION MATRIX 

DEFINE THE ELEMENTS FOR FIRST AND LAST COLUMNS 



NEW086 50 
NEW0866 0 
NEW086 7 0 
NEW08680 
NEW086 90 
NEW087 0 0 
NEW087 1 0 
NEW08720 
NEW08730 
NEW0874 0 
NEW0875 0 
NEW0876 0 
NEW0877 0 
NEW08780 
NEW08790 
NEW088 0 0 
NEW088 1 0 
N EW08820 
NEW08830 
NEW088 A 0 
NEW08850 
NEW0886 0 
NEW0887 0 
NEW0888 0 
NEW0889 0 
NEW089 0 0 
NEW08910 
N EH08920 
NEW08930 
NEW08940 
NEW08950 
NEW0896 0 
NEW08970 
NEW08980 
NEW08990 
NEW09000 
NEW09010 
NEW09020 
NEW09030 
NEW09040 
NEH09050 
NEW0906 0 
NEW0907 0 
NEW09080 
NEW09090 
NEW09100 
NEW09110 
NEW09120 
NEW09130 
NEW09190 
NEW09150 
NEW09160 
NEW09170 
NEW09180 
NEW09190 
NEW0920 0 
NEW0921 0 
NEW09220 
NEW09230 
NEH0924 0 
NEW09250 
NEW0926 0 
NEW09270 
NEW09280 
NEW09290 
NEW09300 
NEW0931 0 
NEW09320 
NEW09330 
NEN0934 0 
NEW09350 
NEW0936 0 
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DO 170 1=1,10 
DO 171 J = 0,I 

LPHPTR(I,J,0)=(((-1.0)*ALPHA)#XJ) 

IMJ=I-J 

LPHPTRd,lMJ,I) = LPHPTRd,J,0) 

171 CONTINUE 

170 CONTINUE 

C FIRST ORDER MATRIX HAS DFEN DEFINED ABOVE 

C SECOND TO TENTH ORDER BELOW 
DO 175 K = 2, 1 0 
KM1 =K-1 

DO 176 L = 0 , KM1 
DO 177 M=1,KM1 
MM1=M-1 

LPHPTR(K,L,M)=LPHPTR(KMl,L,M)+(( L PHPTRCKMl , L , MM1 ) ) *(-l .0D+00) 
CXALPHA) 

177 CONTINUE 

176 CONTINUE 

DO 178 N = 1 , KM1 
KMN=K-N 

LPHPTR(K,K,N)=LPHPTR(K,0, KMN) 

178 CONTINUE 
175 CONTINUE 

C xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

C CALCULATE THE FINAL FILTER 

C FOR HIGHPASS, THE PROTOTYPE MUST SUBSTITUTE (-Z) FOR (Z) 

IFCPASS. EQ.2)THEN 
DO 180 K = l,2 
DO 181 L = 0, ORDER 

DPROTCK, L)=DPROT(K, L)*( (-1 . 0)**L ) 

181 CONTINUE 

180 CONTINUE 

ENDIF 

DO 185 K = l,2 
DO 186 L=0, ORDER 
SUMR=0 . 0D+00 
DO 187 M=0, ORDER 

TEMP=DPROT( K,M)*LPHPTR( ORDER, M,L) 

SUMR=SUMR+TEMP 
187 CONTINUE 

DFLTRC K, L ) = SUMR 
186 CONTINUE 

185 CONTINUE 
ENDIF 

C XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C BANDPASS AND BANDSTOP FORMULAS 

I F ( ( PASS . EQ . 3 ) . OR . ( PASS . EQ . A ) )THEN 
C LOAD THE BANDPASS/ BANDSTOP TRANSFORMATION MATRIX 

C FIRST LOAD THE FIRST AND LAST COLUMNS 

DO 225 1=1,10 
DO 226 J=0,I 
BPBSTRd, J, 0)=C*XJ 
12=2*1 
IMJ =1 - J 

BPBSTR(I,IMJ,I2)=BPBSTR(I,J,0) 

226 CONTINUE 
225 CONTINUE 

C LOAD THE REST OF THE FIRST ORDER MATRIX 
BPBSTRd, 0, 1 )=(-l . 0D+00)*B 
BPBSTRd, 1,1) = (-1 . 0D+0 0)*B 

C LOAD THE 2ND-10TH ORDER TRANSFORMATION MATRICES 
DO 230 K=2, 1 0 
K2=2*K 
K2M1=K2-1 
KM1 = K-1 
K2M2=K2-2 
KM1 2 = 2*KM1 
KM1 2M1 =KM1 2-1 

C LOAD THE SECOND AND SECOND TO LAST COLUMNS 
DO 231 L =0 , KM1 

BPBSTRCK,L,1)=BPBSTR<KM1,L, 1)+CBPBSTRCKM1,L,0)X(-1 .OD+00)XB) 
BPBSTR(K,0,K2M1)=(BPBSTR(KM1,0,KM12)X(-1 ,OD+00)XB)+ 



NEW09370 
NEW09380 
NEW0939 0 
NEW09A00 
NEW09A1 0 
NEW09A20 
NEW09A30 
NEW09AA0 
NEW09A50 
NEW09A60 
NEW09A70 
NEW09A80 
NEW09 A9 0 
NEW09500 
NEW09510 
NEW09520 
NEW09530 
NEW095A0 
NEW09550 
NEW09560 
MEMO 9 57 0 
NEW09580 
NEMO 9 590 
NEW09600 
NEW096 1 0 
NEW09620 
NEMO 96 30 
NEW096A0 
NEMO 96 50 
NEW09660 
NEW0967 0 
NEW0968 0 
NEW 096 90 
NEMO 9 7 00 
NEW097 1 0 
MEMO 97 20 
NEW09730 
NEMO 97 AO 
NEW09750 
NEW0976 0 
NEW 09 77 0 
NEMO 97 80 
NEW09790 
NEMO 9 8 00 
NEW09810 
NEW09820 
NEW 098 30 
NEW 098 AO 
N EW 098 50 
NEW09860 
NEW09870 
NEW0988 0 
NEW09890 
NEMO 99 0 0 
NEMO 99 10 
NEMO 9 9 20 
NEW09930 
NEW099A0 
NEW 09 950 
NEW0 996 0 
NEW09970 
NEW09980 
NEW09990 
NEW10000 
NEW10010 
NEW10020 
NEW10030 
NEW100A0 
NEW10050 
NEW10060 
NEW10070 
NEW10080 
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231 

C 



237 

236 

C 



238 

230 

C 

C 

C 



241 

240 



247 

246 

245 

C 

C 

C 



C 



251 

250 



300 



301 



C 



C (BPBSTR(KM1,0,KM12M1)*C) 

KML=K-L 

BPBSTRCK,KML,K2M1)=BPBSTRCK,L,1) 

CONTINUE 

NOW LOAD INTERIOR ELEMENTS EXCEPT FOR LAST ROW 
DO 236 L = 0 / KM1 
DO 237 M=2,K2M2 
MM1=M-1 
MM2=M-2 

BPBSTR(K,L,M)=BPBSTR<KM1,L,M)+<BPBSTR(KM1,L,MM1)X<-1 .0D+00)X 
C Bl+CBPBSTRCKMl, L,MM2)XC) 

CONTINUE 

CONTINUE 

NOW FILL THE LAST ROW 
DO 238 N=1,K2M1 
K2MN=K2-N 

BPBSTR(K,K/K2MN)=BPBSTR(K/0,N) 

CONTINUE 

CONTINUE 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

CALCULATE THE FINAL FIlTFRv 

FOR THE BANDPASS CASE REPLACE (2) WITH (-Z) 

I F C PASS . EQ. 3 ) THEN 
DO 240 M=l,2 
DO 241 N-0 , ORDER 

DPROTCM,N)=DPROTCM,N)XCC-l .0D+00)XXN) 

CONTINUE 
CONTINUE 
ENDIF • 

ORDER2=2*ORDER 
DO 245 K = 1 , 2 
DO 246 L = 0 , 0RDER2 
SUMR = 0 . 0D+ 00 
DO 247 M=0, ORDER 

TEMP = DPROT C K , M ) XBPBSTRC ORDER, M, L ) 

SUMR=SUMR+TEMP 
CONTINUE 
DFLTRC K, L ) =SUMR 
CONTINUE 
CONTINUE 
ENDIF 

xxxxxxxxxxxxxxxx 
OUTPUT OPTIONS 
xxxxxxxxxxxxxxxx 

I F( ( PASS . EQ . 3 ) . OR . ( PASS . EQ . 4 ) )THEN 
Ell D=2 BORDER 
ENDIF 

IF( (PASS. EQ. 1) .OR. (PASS. EQ.2)) THEN 
END=ORDER 
ENDIF 

MAKE DENOMINATOR COEFFICIENT OF Z**END EQUAL TO 1 
DO 250 M=l,2 
DO 251 N = 0 , END 

DFLTRC M, N) = DFLTRC M, N)/DFLTRC 2 , END) 

CONTINUE 

CONTINUE 

CALL EXCMS ( ’CLRSCRN* ) 

WRITEC 6 , X) 

MRITEC6 , XJ * THE FINAL DIGITAL FILTER IS: 1 
WRITEC 6, x) 

WRITEC 6 , X) * NUMERATOR COEFFICIENTS: • 

DO 300 I = END, 0,-1 
WRITEC 6 , 902 ) DFLTR(1,I),I 
CONTINUE 

WRITEC 6 , X ) • DENOMINATOR COEFFICIENTS: 1 
DO 301 J=END,0,-1 
WRITEC6/902) DFLTRC2,J),J 
CONTINUE 
WRITEC 6 , X ) 

WRIT E C 6 / x ) ' xx ENTER 1 TO CONTINUE XX* 

READC5,X) LOOK 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



NEW10090 
NEW10100 
NEW10110 
NEW10120 
NEW10130 
NEW10140 
NEW10150 
NEW10160 
NEW10170 
NEW10180 
NEW10190 
NEW10200 
NEW10210 
NEW10220 
NEW10230 
NEW10240 
NEW10250 
NEW10260 
NEW10270 
NEW10280 
NEW10290 
NEW10300 
NEW10310 
NEW10320 
NEW10330 
NEW10340 
NEW10350 
NEW10360 
NEW10370 
NEW10380 
NEW10390 
NEW104Q0 
NEW10410 
NEW10420 
NEW10430 
NEW10440 
NEW10450 
NEW10460 
NEW10470 
NEW10480 
NEW10490 
NEW10500 
NEW10510 
NEW10520 
NEW10530 
NEW10540 
NEN10550 
NEW10560 
N EMI 0 57 0 
NEW10580 
NEW10590 
HEW10600 
N EMI 06 10 
NEM10620 
NEW10630 
NEW10640 
NEW10650 
NEM10660 
NEW10670 
NEM10680 
N EMI 06 9 0 
N EMI 07 00 
NEW10710 
NEW10720 
NEW10730 
NEW10740 
NEW10750 
NEW10760 
NEW10770 
NEW10780 
NEW10790 
NEW10800 
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C 



310 



320 



OPTION TO SPOT CHECK FILTER OUTPUT MAGNITUDE AND FREQUENCY 
CALL EXCMS C'CLRSCRN') 

WRITEC6,*) 'THE TRANSFER FUNCTION FOR THIS FILTER HAS BEEN' 
WRITEC6,X) 'COMPUTED. DO YOU WANT TO SPOT CHECK THE MAGNITUDE' 
WRITEC6,*) 'AND PHASE OF THE FILTER OUTPUT FOR A GIVEN * 
WRITEC6,X) 'FREQUENCY?' 

WRITEC6,X) ' XX I-YES/2-NO XX • 

READ(5,X) LOOK 
I F( LOOK . EQ . 1 )THEN 
CALL EXCMS ( 'CLRSCRN ' ) 

I F( AFORDF . EQ . 1 )THEN 

WRITEC6,X) 'WHAT ANALOG FREQUENCYCIN KHZ)' 

WRITE(6,X) ' DO YOU WANT TO CHECK?' 

WRITE(6,X) • XX DECIMAL XX' 

READ(5,X) F 
DF=(F/SSFREQ)X2.0XSPI 
ENDIF 

IFCAF0RDF.EQ.2)THEN 

WRITEC6/X) 'WHAT DIGITAL FREQUENCY DO YOU WANT TO CHECK?' 
WRITE(6,X) • XX DECIMAL XX » 

READ(5,X) DF 
ENDIF 
X=COS ( DF ) 

Y=SIN(DF) 

Z=CMPLX( X/ Y) 

XN=REALCDFLTRC1,0)) 

YN = 0 . 0 



XD=REALCDFLTRC2,0)) 

YD= 0 . 0 

HOFZN=CMPLXCXN,YN) 

HOFZD=CMPLXCXD,YD) 

DO 320 N=1 , END 
ZN=REAL(DFLTR(1/N) )X(ZxxN) 
ZD=REAL(DFLTR(2,N))X(2XXN) 

HOFZN=HOFZN+ZN 

HOFZD=HOFZD+ZD 

CONTINUE 

HOFZ = HOFZN/HOFZD 

FMAG=SQRT((REAL(H0FZ)XX2)+(AIMAG(H0FZ)XX2)) 

FMAGDB=20 . OXAL OG1 0 ( FMAG) 

FPHSE=ATAN(AlMAG(HOFZ)/REAL (HOFZ) ) 

MRITEC6 , x) 'THE FILTER OUTPUT WILL BE:' 

WRITEC6,91A) FMAGDB 
WRITEC6,915) FPHSE 
NR I T E ( 6 / x ) 

NR I TEC 6 / X) 'DO YOU WANT TO CHECK ANOTHER FREQUENCY?' 
WRIT EC 6 , x) ' xx 1-YES/2-N0 XX ' 

READC 5 , X) LOOK 
IFCLOOK. EQ.DGO TO 310 
ENDIF 



xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 



OPTION TO LOAD FREQUENCY RESPONSE INTO A DATA FILE 
PAI I FYPMS ( Tl RSPRJJ * > 

1 DO YOU WANT THE FILTERS FREQUENCY RESPONSE ' 

DATA TO BE ENTERED INTO A DATA FILE FOR USE' 

BY A PLOTTING PROGRAM? ' 

NOTE t THE FILE NAME WILL BE FILTER, FILE TYPE' 

WILL BE DATA, 100 POINTS WILL BE CALCULATED' 

FOR A DIGITAL FREQUENCY FROM 0 TO PI, THE FIRST' 
COLUMN WILL BE THE DIGITAL FREQUENCY (EXPRESSED 
AS A FRACTION OF THE SAMPLING FREQUENCY 0-1/2) ' 
THE SECOND COLUMN WILL BE THE MAGNITUDE, AND' 

THE THIRD COLUMN WILL BE THE PHASEC DEGREES ) ' 

XX 1-MAKE DATA FILE/2-D0N"T BOTHER XX* 



WRITEC6 , X) 

WRITEC6, x) 

WRIT EC 6 , X ) 

WRITEC6 , x) 

WRITEC6, x) 

WRITEC6 , x) 

WRITEC6 , X) 

WRITEC 6 , x ) 

WRITEC6,X) 

WRITEC 6 , x ) 

HRITEC6,x) 

READC 5, X) LOOK 
IFCLOOK. EQ.DTHEN 
DO 350 1=0,100 
DF=REALCI)XC5. 0/1000, 
DFR=REALCI)XCSPI/100, 
X = COS C DFR ) 

Y=SINC DFR) 

Z=CMPLXCX, Y) 



0 ) 

0) 



NEW10810 
NEW10820 
NEW10830 
NEW108A0 
NEW10850 
NEW10860 
NEW10870 
NEW10880 
NEW10890 
NEH10900 
NEW10910 
NEW10920 
NEW10930 
NEW109A0 
NEW10950 
NEW10960 
NEW10970 
NEW10980 
NEW10990 
NEW11000 
NEH11010 
NEW11020 
NEW11030 
NEW11040 
NEW11050 
NEW11060 
NEW11070 
NEW11080 
NEW11090 
NEW11100 
NEW11110 
NEW11120 
NEW11130 
NEW111A0 
NEW11150 
NEW11160 
NEW11170 
NEH11180 
NEW11190 
NEW11200 
NEW11210 
NEW11220 
NEW11230 
NEW112A0 
NEW11250 
f J EH 1126 0 
NEW11270 
NEN11 280 
NEW11290 
N EMI 1300 
NEW11310 
NEW11320 
NEW11330 
NEW113A0 
NEW11350 
NEW11360 
NEW11370 
NEW11380 
NEW11390 
' NEW1 1 A 0 0 
NEW11A10 
NEW11A20 
NEW11A30 
. NEW11440 
NEW 1 1 AS 0 
NEW11A60 
NEW11A70 
NEW11A80 
NEW11A90 
NEW11500 
NEW11510 
NEW11520 
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FILE* DFCADD FORTRAN A 



351 



350 

C 

901 

902 

903 
909 
905 
919 

915 

916 



XN = REAL C DFLTRC 1 , 0 ) ) 

YN = 0 . 0 

XD=REAL ( DFLTRC 2/ 0 ) ) 

YD = 0 . 0 

HOFZN=CMPLX(XN,YN) 

HOFZD=CMPLXCXD,YD) 

DO 351 N=1 , END 
ZN=REALCDFLTRC1,N))*CZX*N) 

ZD=REALCDFLTRC2,N))*(Z**N) 

HOFZN=HOFZN+ZN 

HOFZD^HOFZD+ZD 

CONTINUE 

HOFZ=HOFZN/HCFZD 

FMAG=SQRT C C R EAL C HOFZ )x*2 ) + ( AIMAGC HOFZ) **2 ) ) 

IFC REAL ( HOFZ) . NE . 0 • 0)THEN 
FPHSER=ATANC AIMAGC HOFZ)/REALC HOFZ)) 

ELSEIFC REAL C HOFZ) .EQ.O .0)THEN 
I FC AIMAGC HOFZ) . GT . 0.0 )FPHSER=1 .570796 
I FC AIMAGC HOFZ) .LT. 0.0) FPHSER=-1 .570796 
IFC AIMAGC HOFZ) .EQ.O . 0 ) FPHSER=0 . 0 
ENDIF 

FPHSE=FPHSERXC180.0/SPI) 

MRITEC 8,916) DF, FMAG, FPHSE 
CONTINUE 
ENDIF 

XXXXX*XX**XXX*X*X**KXXXXXXXX**X**XX*X*XKKXXXK***********X*X***** 

FORMATC 5X, FI 0 . 5, ' S**',I1) 

FORMAT C 5X> FI 0 . 5, 1 Z**',I1) 

FORMATC3X, *ALPHA= »,F9.9) 

FORMATC 3X, ' K= *,F9.9) 

FORMATC3X, 'THE LOMPASS DIGITAL PROTOTYPE FREQ. IS »,F5.3) 
F0RMATC3X, 'MAGNITUDE = »,F8.2,' DB ' ) 

FORMAT C 3X, ' PHASE = ',F6.2,' RADIANS') 

FORMAT (3X,F9.3,10X,F5.3,10X,F6.2) 

STOP 

END 



NEM11530 
NEM11590 
NEM11550 
N EMI 1560 
NEW11570 
NEM11580 
NEH11590 
N EMI 1 6 0 0 
NEM11610 
NEM11620 
NEM11630 
NEM11690 
NEM11650 
NEM11660 
NEM11670 
NEM11680 
NEM11690 
NEM11700 
NEM11710 
NEM11720 
NEM11730 
NEM11790 
NEM11750 
NEM11760 
NEM11770 
N EM 1 1 7 8 0 
NEM11790 
NEM11800 
N EM 1 1 8 1 0 
NEM11820 
NEM11830 
NEM11890 
NEM11850 
NEM11860 
NEM11870 
NEM11880 
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